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Chapter 1 

Introduction 


The following report describes the research activity performed by the Associate during the 
period starting January 14th, 2002 until August 16th, 2004 at the NASA John C. Stennis 
Space Center (SSC) in Mississippi. The Associate wishes to thank his advisor, Dr. Fernando 
Figueroa; Dr. Enrique Barbieri, a summer visiting faculty from Tulane University; Dr. Bill 
St.Cyr, Dr. Bob Field and Mr. Jared Sass, NASA engineers; and Ms. Jamie Granger, a summer 
undergraduate student from Tulane University for their valuable cooperation in the project. 

1.1 Objectives, Scope and Methodology 

A rocket test stand and associated subsystems are complex devices whose operation requires 
that certain preparatory calculations be carried out before a test. In addition, real-time control 
calculations must be performed during the test, and further calculations are carried out after a 
test is completed. The latter may be required in order to evaluate if a particular test conformed 
to specifications. These calculations are used to set valve positions, pressure setpoints, control 
gains and other operating parameters so that a desired system behavior is obtained and the 
test can be successfully carried out. Currently, calculations are made in an ad-hoc fashion and 
involve trial-and-error procedures that may involve activating the system with the sole purpose 
of finding the correct parameter settings. The goals of this project are to develop mathematical 
models, control methodologies and associated simulation environments to provide a systematic 
and comprehensive prediction and real-time control capability. The models and controller 
designs are expected to be useful in two respects: 

1. As a design tool, a model is the only way to determine the effects of design choices without 
building a prototype, which is, in the context of rocket test stands, impracticable. 

2. As a prediction and tuning tool, a good model allows to set system parameters off-line, so 
that the expected system response conforms to specifications. This includes the setting 
of physical parameters, such as valve positions, and the configuration and tuning of any 
feedback controllers in the loop. 

This work fits into a framework of current efforts in developing a set of tools that simplify 
prediction tasks routinely done at the test stands. The tools are user-friendly software that 
incorporates knowledge about particular systems and uses it to carry out predictions of system 
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operation under diverse conditions. Building such knowledge is the most important part of the 
work and involves mathematical modeling, data collection, model validation and analysis of 
model properties. Work on a graphical user interface (GUI) for the E-l hydrogen mixer is also 
considered a project goal. 

The following specific activities have been carried out as part of the research project: 

1. Literature review and consultation with NASA engineers. 

2. Construction of mathematical models. 

3. Data collection and organization. 

4. Use and development of simulation programs. 

5. Model validation with experimental data. 

6. Dynamic model analysis and controller synthesis. 

7. Closed-loop control simulations and feasibility studies. 

8. Work in control-theoretical problems inspired by practical problems. 

9. Preparation of technical articles and presentations. 

10. Programming of specialized simulation environments. 

The work methodology roughly follows the sequence of activities presented above. Some itera- 
tions are often required for model and control law refinement. 

1.2 Summary of Accomplishments 

The tenure period has been very productive, as evidenced by the number of technical papers 
resulting from the research project. The following list of achievements summarizes the main 
points: 

• Construction and experimental validation of a model for the E-l Hydrogen Mixer. 

• Dynamic analysis: equilibrium, controllability, linearizability. 

• Simulation tools in Matlab/Simulink. 

• Controller designs: Small-signal and feedback linearization, Sliding Mode. 

• Theoretical work on constrained Sliding Mode Control (new contribution to the field). 

• Design and programming of Matlab-based Graphical User Interface. 

• Preliminary studies in Run Tank dynamics. 

Seven conference presentations [3, 21, 22, 25, 26, 27] -and corresponding proceedings articles- 
have resulted from the project, along with one accepted journal article [23] and one in the 
submission process. In addition, the work in the Graphical User Interface has resulted in a 
NASA Software Disclosure, to be reported in the NASA Tech Briefs [24] . 
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Chapter 2 

The E-l Hydrogen Mixer 


The hydrogen mixer subsystem of the E-l test facility is used to condition the flow of liquid 
hydrogen (LH2) before it is fed into a test article. Specifically, the temperature of the liquid 
at the exit of the mixer must be controlled, along with the pressure inside the mixer and the 
total output flow rate. This is accomplished in the mixer by injection of a stream of gaseous 
hydrogen (GH2). The rate at which GH2 is injected determines both the transient and steady- 
state temperature behavior of the conditioned propellant. The system is depicted in Fig. 2.1 
below. 

2.1 Assumptions - Preliminary Model 

A preliminary model is derived based on the following assumptions: 

• Thermodynamic properties of LH2 and GH2 at the mixer inputs are constant. 

• The mass flow rates of LH2, GH2 and output mixture are independent control variables 
that can be manipulated at will. 

• The mixer operates adiabatically (i.e. , is thermally insulated from the surroundings) and 
isobarically (the pressure drop in the mixer is negligible). 

• Thermodynamic properties of the exit stream (i.e., temperature) equal those inside the 
mixing chamber. Also, GH2 condenses completely inside the mixer, rejecting heat to the 
LH2 and thus increasing its temperature. There is no two-phase flow at the exit. 

• The output flow is incompressible. 

2.1.1 Governing equations 

The two primary equations used to model the process in the mixer are the conservation of 
mass and the conservation of energy (first law of thermodynamics for a control volume). Both 
equations are used in their more general forms, to account for transient effects. The conservation 
of mass equation is 

rh cv = ^ rhj - ^ m e (2.1) 
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Figure 2.1: Mixer Subsystem 


where m cv is the mass inside the mixer (control volume) and rhj and rh e are the input and 
exit mass flows, respectively. In this case, there are two input and one exit mass flows. Under 
transient conditions, the time derivative of m cv is nonzero. Denoting the LH2 input mass flow 
as rhi , the GH2 mass flow as rh g , and the exit mass flow as m e , the mass balance equation 
becomes 

fn cv = fill + da g — ih e (2-2) 

The first law of thermodynamics for a control volume is stated as: 

r ]TT y2 y2 

^ rhi(hi 4 — ^ — h gZi) — fn e (h e 4 — 4- gZ e ) + Q cv — W cv (2.3) 

Here U cv is the (internal) energy of the mass inside the mixer, hi and h e are the specific en- 
thalpies of the input and exit flows, and V t , V e , Zi and Z e denote their velocities and elevations. 
Since the mixer does not introduce important velocity changes, the kinetic energy terms will be 
neglected, and so will be the changes potential energy due to elevation. The rate of heat trans- 
ferred through mixer walls, Q cv , and the mechanical work done inside the mixer are considered 
zero. The derivative of the internal energy of the mixer can be expanded as: 

dU cv d(m cv u cv ) 

77 — 77 — ///,(■;■ //(■;■ T Ucv'Tll’CV 

at at 

where u cv is the specific internal energy of the mixer. Following the assumption of constant 
input properties, the energy equation becomes 

m cv Ucv + u cv m cv = rhihi + m g h g - m e h e (2.4) 

where hi and h g are the constant input enthalpies and h e is the exit enthalpy, which will 
vary under transient conditions. If the three mass flows are specified, Eqs. (2.2) and (2.4) 
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form a system of ordinary differential equations on the variables m cv and u cv . However, the 
exit enthalpy h e is unknown, and therefore a solution cannot be obtained. An additional 
assumption is therefore necessary to completely determine the system. It will be assumed that 
the thermodynamic properties inside and at the exit of the mixer are the same. Specifically, the 
exit internal energy u e is equal to the mixer internal internal energy u cv . Using the definition 
of enthalpy and differentiating: 


dh e = du e + dp e v e + p e dv e 

Using the constant pressure and incompressibilty assumptions gives 

he — U ( — 'dev 

therefore 

h e = u cv + C (2.5) 

where C is a constant determined from the initial conditions. Now, Eqs. (2.2), (2. 4), and (2.5) 
can be solved for any set of initial conditions. Note that the initial values of m cv and u cv are 
not easily measurable in the real process, but that can be indirectly calculated from the input 
and output properties. 

2.1.2 The mixer as a control system 

The mixer equations must be put in a form adequate for control system analysis and design. 
One of such forms is the state-space representation, which has a convenient mathematical form. 
Let the independent control inputs be defined as: 

Vh = mi 
i>2 = dig 

= m e 

Let the two states be denoted as: 

X\ — 'Ucv 
X2 — TYlcv 

Using the governing equations and the above definitions, the mixer control system can be 
represented as: 

xi = —hii/)i + hgi/)2-Ci/>3-xi(il)i + il>2) ( 2 . 6 ) 

X'2 

X 2 = V’i + Y>2 - ^3 (2.7) 

In vector form, the mathematical model falls in the category of a 2-state, 3-input bi-linear 
model of the form 


X\ = f T (x 1 ,x 2 )'tp 
x 2 = g (V>) 
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where ijt = [ipi 1 P 2 ' l p 3 ] T is the vector of control flows and the functions / and g are given by 


f(x i,x 2 ) = 


hi-x i 
X2 

hg~X 1 

X2 

_C_ 

X2 


V>2, fo) = Ipl + ^2 - V>3 


(2.8) 


(2.9) 


Preliminary observations indicate that the methods of feedback linearization [6], backstepping 
control [13, 29] or differential flatness [18] could be candidates for controller designs. 

2.1.3 Control valve models 

The following expression models the volumetric flow of liquid through a valve: 


Q = C 


CvV^n Pout 

TIT 


where C is a constant dependent on choice of units, C v is the valve flow coefficient, S w is 
the specific gravity of the liquid relative to water, and P{ n and P ou t are the inlet and exit 
pressures, respectively. The specific gravity is taken constant, at a reference temperature. For 
compressible fluid, the following formulas are considered: 


Q = { 


CC V \! Pi, i^ out ,when P in < 2 P out 


C‘ 


7 CyP in 
P°RS a 


,when Pj n > 2 P t 


out 


Here, °R is the fluid temperature in degrees Rankine, S a is the specific gravity of the gas relative 
to the air, taken constant at 14.7 psia and 60°-F and C, C’ are constants depending on the choice 
of units. Note that the flow given by these expressions is in Standard Cubic Feet per Second 
(SCFS), that is at 14.7 psia and 520°i7. A modification must be introduced to adjust for actual 
flow conditions using the ideal gas law. The actual flow Q a is computed as 

= JLq 

^ P 520 V 

whereas the mass flow is obtained by multiplying Q a by the density of the fluid matching P and 
T. The P, T , p combination may be taken at valve inlet or outlet, the result being theoretically 
the same. For the case of hydrogen, the formulas for mass flow become: 


For LH2 flow: 


IT = 1.76 x 10 ~ 2 C v y/{P m - P out )p 


(2.10) 
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For GH2 flow: 


2.857 x 10~ 2 C V J P 2 n - P 2 ut ,when P in < 2 P c 


out 


w = < 


2.423 x 10 ~ 2 C v P in ^ 


,when Pi n >2 P c 


out 


( 2 . 11 ) 


In the above formulas, the mass flow is given in lb/ s, the pressures in psi , the densities in lb/ ft 3 , 
and the temperature in °R. The value used for the specific gravity of hydrogen is Sa = 0.0699 
(relative to air at 1 atm and 70°.F) and the density of water in Eq. (2.10) has be taken as 
62.4 Ibm/ ft 3 . 


2.1.4 Example 

The following numerical example compares the values given by the above formulas with the 
ones in [2], where the formula for LH2 was used for GH2. Let Pi n = 13500 psia, P ou t = 6000 
psia, Ti n = 90°F, T out = 135 °F, pi n = 2.9 Ibm/ ft 3 and p ou t = 1.53 Ibm/ ft 2 . Let the flow 
coefficient for a fully open valve be 230 and let the valve be 4.35 percent open. This gives 
C v = 10. Since > 2 P ou t, the second of Eq. (2.11) is used, yielding a mass flow of 16.54 
lbrn/s when inlet conditions are used, and of 20.36 lbrn/s when oulet conditions are used. If the 
formula for liquid is used, the mass flow is 26 lbrn/s. It should be noted that there are more 
sophisticated gas flow formulas which employ empirical correlations from the manufacturer. At 
this point, however, the above formulas will be considered to provide sufficient information for 
control analysis and design purposes. 


2.2 Revised mixer model 


A more realistic model can be derived by relaxing the constant liquid density and mixer pressure 
assumptions. This will require, however, the inclusion of a functional thermodynamic relation, 
which increases the complexity of the model. The mixer pressure as a function of density 
and temperature will appear either as a table look-up, or explicitly as a correlation formula. 
Modeling assumptions are the same as in the preliminary model, except that the density at 
the output and mixer pressure are allowed to vary. Also, the control flows are related to fluid 
properties and valve characteristics through valve models. The control variables now become 
the valve C v coefficients, which in turn can be ultimately related to stem position and voltage. 
The latter is expected to be the controlled variable in an actual system. Considering that 
the output flow is entirely in the liquid phase and using the above valve models, the mass 
conservation equation becomes: 


Pcv ~ V 


CiC vl y/(Pi - P)pi + C 2 C ^Y ] VT 9 P 9 ~ C 3 C vey /(P - P s )p c 

-Lfl 


(2.12) 


where V is the constant mixer volume, Ci, C 2 and C 3 are constants, C v i, C vg and C ve are the 
valve flow coefficients for the LH2, GH2 and exit flow control valves, Pp P g are the pressures of 
the LH2 and GH2 upstream of the input control valves, and P is the mixer pressure, assumed 
to be varying in time but uniform throughout the mixer. The function f(P g ,P) is chosen 
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according to Eq. (2.11). Note that it has been assumed that the exit density equals the density 
of the mixer. The energy conservation equation becomes: 


'U'cv — 


Vpc 


CiC vl ^(P l - P)p l h l + C 2 C V9 f ^ P) VT 9P9 h 9 - C 3 c ve y/(P - p s )p cv h cv - 


Tflcv 'U' c 


(2 - 13) 

Again, the exit enthalpy and density are considered equal to those within the control volume. 
The enthalpy at the exit can be related to other variables: 


, , P P 

hcv — h e — u e T — Ucv T 

Pe Pcv 


Also, the mixer pressure is a function of the density and internal energy: 


P — P{^CVl Pcv) 


This function can be implemented as a table look-up or as an empirical correlation. For simu- 
lation purposes, the time variation of Pi, P g , P s , hi, h g and T are supposed to be known. As 
the model is expanded to include the run tank dynamics, some of these variables will be put 
in terms of quantities that are measurable in the real system. For analysis purposes, the mixer 
model is a system of two ordinary differential equations with three independent inputs (the 
valve flow coefficients). The equations become: 


Pcv 


1 

V 


C 1 C vl ^/(P l ~P(u cv ,p cv ))p l + C 2 a g /(Fg,P( p CT,Pci,)) sjT gPg - C 3 C ve ^/(P(u cv ,p cv )-P s )p c 

Ffi 


'U'cv — 


Vp c 


C^C^Pi-Piu^Pa^piihi - U CV ) + C 2 C vg /(Pg ’ P( p CT ’ PCT)) y/T g p g [h g - U CV ) 

* a 


n 1 ocr 1 r* P( u cvi Pcv) Ps p/ \ 

0.\S^>C^C ve \ P(u cv ,Pcv) 


(2.14) 


2.3 Mixer Model: Static Properties 

For controller design purposes, it is of interest to gain insight on the mathematical properties of 
the model. Model structure can be often exploited to design better controllers. In particular, 
the number and nature of the equilibrium points should be investigated, as well as dynamic 
properties such as linearizability, stability of internal dynamics, etc. 


2.3.1 Characterization of the Equilibrium Point 

The equilibrium points are defined as the constant values of the state which result in zero 
derivatives for constant control inputs. In this case, the control inputs are the C v coefficients. 
The equilibrium point indicates the steady values of the density and internal energy of the mixer 
when the control valves are set at fixed positions, for a given set of input flow characteristics. 
For given values of the input fluid properties and C v coefficients, setting p cv = 0 establishes 
that m e = m g + rip and results in an expression relating the equilibrium density p cv to the 
equilibrium mixer pressure P: 

Pcv = p(P) (2.15) 
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Graphical Interpretation of Mixer Equilibrium 



Density, kg/m 3 

Pressure, MPa 


Figure 2.2: Mixer equilibrium 


Setting u cv = 0 gives the equilibrium exit enthalpy h e = h cv in terms of the input enthalpies 
and the mass flows: 

- mihi + m a h 0 . . 

h cv = , . 9 9 = h{P) 2.16 

mi + m g 

The above enthalpy must match the thermodynamic property data at P and p cv of Eq. 2.15, 
that is 

hcu — hfh ( Pcv j -P) 

Substituting Eqs. 2.15 and 2.16 into the above equation results in a single expression which 
gives the equilibrium pressure: 

h(P) = h th (p(P),P) (2.17) 

A graphical interpretation of the equilibrium solution is shown in Fig. 2.2. The curve pP 
vs. P has been drawn on the base plane. This plane curve is mapped to a space curve by 
the thermodynamic property function hfj-, . The equilibrium exit enthalpy h(P) is a function 
of pressure only and the corresponding surface has also been graphed. The point where the 
space curve pierces the surface is the equilibrium point of the system. It is seen to happen 
approximately at P = 5800 psia and p cv = 3.5lbm/ ft 3 . The monotonicity of the curve and 
surface suggests that there exists only one equilibrium point in the general case. 

2.3.2 Iterative solution for the equilibrium point 

An iteration scheme can be used to find the exact value of the equilibrium, perhaps using a 
graph to first find an approximate solution, as above. One possible iteration method is given 
below: 
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1. Start with an initial pressure guess, Pq. 

2. Calculate the mixer equilibrium density p cv from Eq. (2.15). 

3. With P 0 and p cv calculate the enthalphy h t h using thermodynamic data. 

4. With P 0 calculate the steady-state enthalpy h cv from Eq. (2.16). 

5. Calculate the error e = h cv — hth- 

6. Correct the pressure guess using P <— P + ke, where A; is a convergence factor. 

7. Repeat steps 2 through 6 until the error is within a prescribed tolerance. 

The above steps have been coded in Matlab and it has been determined that, for the expected 
values of mixer operation, k = 0.1 is a good convergence factor. The following is a sample run 
using the pressure guess value from the graphical method: 

» solv_mxr_eq 
ans = 

Convergence achieved in 23 iterations 
P = 

5 . 8745e+003 


e = 


9 . 6981e-004 


rho_cv = 
3.2538 


h_th = 

609.8934 

2.3.3 Methods of evaluating P(p,u) 

The NIST thermodynamic data for Hydrogen can be accessed using a the MIPROPS program. 
For Matlab simulation purposes, it would be desirable to have the pressure information as a 
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table, since the table look-up blocks included in Simulink greatly simplify the simulation setup. 
However, MIPROPS provides density, internal energy, -and other thermodynamic data- for 
regularly spaced temperature points and a specific pressure. Conversion of such data into a 
rectangular array with pressure entries for given density and internal energy is not a trivial task. 

If pressure and temperature ranges are set, and density and energy data are extracted using 
MIPROPS, it is seen that the set of all density-energy pairs has a non-rect angular boundary 
upon plotting. Moreover, the density-energy points are irregularly spaced. Extending such 
boundary to a rectangular one and regularizing the grid of points to conform a rectangular 
array is cumbersome and introduces error due to interpolation. Instead of attempting to obtain 
a pressure table, it is more convenient to write a function that determines the pressure based 
on available data. Three methods are considered and explained as follows. 

2.3.4 Method 1 

This method considers that the internal energy is a function of temperature only. In reality, 
the energy-temperature curve varies with pressure. The given internal energy is mapped into 
temperature using the curve for some pressure within the expected range of operation. Once 
the temperature has been obtained, a pair of bracketing temperatures are found in the density- 
energy table, and interpolation is used to find the pressure. 

2.3.5 Method 2 

This method finds bracketing energies among all entries of the density-energy table. From the 
bracketing energies found, only those whose corresponding densities bracket the given density 
are kept. Those energy brackets are interpolated for density using the given energy, and the 
bracket that yields the interpolated density that is closest to the given density is selected. The 
pressure corresponding to the bracket is returned by the function. 

2.3.6 Method 3 

This method uses Method 1 to find the pressure once the temperature is known. The tempera- 
ture, however, is not found exclusively from the energy, but from the full density-energy data. 

The above methods have been programmed in the Matlab function files getpress.m, get- 
temp. m and getpt.m listed at the appendix. 

2.3.7 Comparison of Methods 

The Simulink diagrams of the mixer model using Methods 1,2 and 3 are shown in Figs. 2.3, 2. 4, and 2.5. 

A comparison of the methods is shown in Fig. 2.6. It is seen that methods 2 and 3 provide 
results that are closer together relative to method 1. It is also seen that the basic model of 
Section 2.1.1 introduces significant error in energy. However, upon translation of the energy 
into temperature, the error is reduced. It must also be pointed that method 2 has the largest 
execution time, while method 1 and the basic model have the smallest. Method 3 provides 
a good compromise between speed of execution and accuracy. Fig. 2.7 shows how the model 
in conjunction with method 3 can be used to simulate a change in valve position. Initially, 
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Figure 2.3: Mixer Simulation - Method 1 
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Figure 2.4: Mixer Simulation - Method 2 
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Figure 2.6: Mixer Simulation - Comparison of Methods 
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Figure 2.7: Mixer Simulation - Varying Cv 


only LH2 flows through the mixer. Then, the GH2 valve is suddenly opened, resulting in a 
temperature increase of the exit flow. 


2.3.8 Valve positions for a prescribed static operating condition 


Having three independent controls in a two-state model allows the selection of steady values for 
the states, and, in addition, an extra degree of freedom is available. This degree of freedom can 
be used to fix the exit mass flow with the desired thermodynamic properties, as shown next. 
Supppose it is desired to have a given exit mass flow, with prescribed temperature (measured 
at the outlet of the exit valve). Let these quantities be denoted by W e and T. The back 
pressure P s at the outlet valve exit is assumed to be constant. The problem is to determine 
the valve coefficients that achieve this. Two degrees of freedom are used for W e and T, and 
the third one is used to meet a desired mixer operating pressure, P. The back pressure P s and 
T determine the enthalpy h, which is constant across the valve. Therefore, the enthalpy and 
mixer pressure at the exit valve inlet are known and, in turn, determine the mixer density from 
thermodynamic data. A static energy and material balance gives the required input flows that 
achieve a prescribed flow-exit enthalpy combination. The following formulas are straightforward 
to derive: 


Wi 


W 0 


W e {h - hg) 

hi - hg 

Wejhi - h) 

hi - hg 


(2.18) 


where h g and hi are the gas and liquid supply enthalpies, respectively. Using the exit flow 
and the input equilibrium flows from Eq. 2.18, the three valve coefficients are determined from 
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Eq. (2.10) and (2.11). A Matlab program mxr-cv.m carries out the necessary computations. 
As an example, an operating condition is targeted. Then, the resulting C v values are fed into the 
inverse program which computes the operating point. The operating conditions are recovered 
with good accuracy. Suppose the C v values that achieve an exit flow of 80 lb /s at a temperature 
of — 300°.F with a mixer operating pressure of 7000 psia are sought. The values returned are 

» [Cvl , Cvg , Cve] =mxr_cv (-300 , 80 ,7000) 
h = 

384.8930 


Cvl = 


43.6348 


Cvg = 

5.0234 

Cve = 

56.9665 

The equilibrium-finding program returns 
» solv_mxr_eq 

ans = 

Convergence achieved in 111 iterations 
P = 

7 . 0023e+003 


e = 


9 . 8624e-004 
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rho_cv = 


4.3345 


h_th = 
386.4838 


exit_flow = 

80.0122 

Although system controllability will be analyzed rigorously in another section, the ability to 
select the operating point as shown above suggests that the system is indeed controllable. 

2.4 Model Validation 

The mathematical model of the mixer given by Eqs.( 2.14) are valid regardless of the substance 
flowing through it. Preliminary tests of mixer operation at SSC have been performed with 
nitrogen as a working substance. One of such tests is used to verify the prediction capabilities 
of the model, together with its accesory programs giving pressure and temperature from given 
density and internal energy. 

2.4.1 Test Sequence 

The objective of the test is to obtain a desired temperature of the exit flow by means of the 
introduction of GN2 to the LN2 stream in the mixer. In the test configuration, the exit valve 
was replaced by an orifice. This does not represent a difficulty in using the model, since it can be 
assumed that the exit valve in the model remains open at a fixed position. The equivalent exit 
valve C v is found from the experimental data. The mixed flow is released into the atmosphere. 
Valve openings are precalculated so that the mass flows match desired quantities, namely, 
wi = 22 lbm/s, w g = 16 lbm/s and, in steady state, w e = 38 lbm/s. Pressure inside the mixer 
must also be maintained at a desired level, namely P = 600 psig. To achieve this, and as 
starting approach, the pressure is continuosly measured and used in a PID control loop around 
the gas valve only. The test sequence is as follows: 

1. Open the LN2 valve to 4% to chill down pipes. 

2. Pressurize the LN2 source to 1000 psig. 

3. Ramp the LN2 valve to 40%. 

4. Manually command the GN2 valve until a mixer pressure of 600 psig is reached. 

5. Close both valves to 3%. (apparently, the previous steps were done to determine the 
opening of the GN2 valve to obtain 600 psig). 
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6. Ramp the LN2 valve back to 40%. 

7. Ramp the GN2 valve to approx. 40% and close the PID control of mixer pressure. Since 
the valve opening is close to the already determined value, a relatively bumpless transfer 
occurs. 

8. Close the valves and terminate the test. 

The following data where gathered at all times during the test: 

• LN2 and GN2 valve positions. 

• LN2 and GN2 inlet pressures (upstream of the valves). 

• LN2 and GN2 inlet temperatures (upstream of the valves). 

• Mixer pressure. 

• Mixer temperatue. 

Note that the flows were not measured, but calculated from pressure differences, fluid properties 
and valve positions. The duration of the test data is 2400 seconds, with data sampled every 0.02 
seconds. This results in data records containing more than 100000 points. For simulation in a 
PC to be practical, the data must be reduced. Data was processed using Matlab’s resample 
command, which decimates and filters the data according to user specifications. Resampling 
and filtering was applied to yield data records having little more than 1000 points each. 

2.4.2 Model Simulation 

The valve positions and inlet fluid properties of the test are shown in Fig. 2.8. The sudden drop 
in temperature observed for the LN2 was not handled by any of the numerical methods available 
in SIMULINK. Since significant valve activity begins after the drop, where LN2 temperature 
is practically constant, this constant value is used as an input at all times. Note that this is 
possible to do given that the mixer response is in the order of seconds. The other resampled 
data is directly fed to the mixer model using the “From Workspace” blocks in SIMULINK. The 
time horizon for the simulation was restricted to the first 920 seconds, which contain the most 
useful information. It is important to note that some numerical problems were obtained due 
to interaction of SIMULINK’s integration routine ode45 (Runge-Kutta 45) and the property 
program yourprops.m. After several attempts, it was determined that a reliable simulation 
is obtained using the following parameters: 

• Integration Method: Variable-step ODE23 (Bogacki-Shampine). 

• Maximum step size: 0.05 

• Initial step size: 0.001 

• Relative tolerance: 0.001 

• Absolute tolerance: determined by SIMULINK. 

• Refine factor: 2 
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Figure 2.8: Valve positions and inlet fluid properties during validation test. 


2.4.3 Model Calibration 

The data measured during the test shows pressures that fall below zero. This may be at- 
tributable to instrument offset which becomes important when actual pressures are low. It is 
to be noted that the mathematical model will not produce negative absolute pressures. This 
induces a slight discrepancy between measured and simulated values. Also, the current version 
of yourprops.m handles pressures starting at 20 psia. This value was used as atmospheric 
pressure during simulation, and it also constitutes a small discrepancy. The valve models are 
calibrated to reproduce the calculated steady flows given by the test operators. The mixer pres- 
sure of 614.7 psia and temperature of —222 °F corresponds to a density of 17.287 lbrn/s. With 
the flow formula of Eq. (2.10) and the 38 lbm/sec exit flow it is obtained that the equivalent 
C ve of the orifice is 21.2. This was later adjusted to C ve = 19 to match experimental data. 
Similarly, it is obtained that C vg = 1.7 and C v i = 8.35 for 37/ 

2.4.4 Results 

Fig. 2.9 shows the resulting state trajectories and flows, along with an error variable which 
ensures that there were no errors in yourprops.m. Figs. 2.10 and 2.11 constitute the validation. 
In view of the expected discrepancies explained above, it can be said that the agreement is very 
good. 

2.4.5 Immediate Improvements to the Model 

The following improvements can be made to have a better prediction capability (listed in order 
of expected effectiveness): 
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Figure 2.9: Simulated mixer states. 


Simulated vs. Measured Mixer Pressure 



Figure 2.10: Simulated and measured mixer pressure. 
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Simulated vs. Measured Mixer Temperature 



Figure 2.11: Simulated and measured mixer temperature. 


• Expand the pressure range of the yourprops program to include near-zero pressure. This 
will allow to set the atmospheric pressure in the simulation to the correct value. 

• Obtain the true value of the bottom pressure of the mixer that was obtained during the 
test. Shift the data and recalibrate the Cv values. 

• Include first-order lag dynamics in the valve models. 

The model can be further refined by performing flow experiments with the valves in order to 
determine more accurate descriptions. Also, the current model uses valve Cv as a primary input 
variable. Valve dynamics can be further studied to replace the control inputs in the model by 
commanded valve position, which are electrical signals to be used in an actual controller. 
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Chapter 3 


Feedback 


Linearization Analysis 


Of the control techniques that are suitable for the model structure, feedback linearization is the 
simplest and most direct. The technique, however, may have disadvantages. Specifically, the 
complexity introduced by the cancellation of all nonlinear dynamics may be to high for a realistic 
implementation with limited computational resources. The other significant disadvantage is 
that all state measurements are required. Peaking of control signals is also a factor of concern, 
especially when actuator saturation and rate limitations are present. The technique, however, 
is worth examining, since essential features of the model appear during the analysis. 


3.1 Output Definitions 

Mixer operation requires tracking of prescribed output mass flow and temperature profiles. As 
seen in a previous section, the extra degree of freedom of the system can be used to specify the 
mixer operating pressure in addition to temperature and mass rate of the out flow, given that 
the output pressure P s is known and constant. While this is true for equilibrium conditions, it 
will now be assumed that simultaneous tracking of the three quantities is also possible. Note 
that this amounts to specifying a desired output vector whose components are temperature at 
the exit valve outlet, mixer pressure, and exit mass flow rate. All of the above quantities are 
rather complicated functions of the state [p u] T . In order to simplify matters, the following 
observations are made: specifying the temperature at the exit valve outlet determines the 
enthalpy, since P s is known and constant. The enthalpy, in turn, is constant across the valve, 
and therefore matches that of the mixer. If mixer density is specified according to the enthalpy 
in a way such that the resulting pressure matches the desired pressure profile, then density 
can be used as a desired output. Since the prescribed density and enthalpy also determine the 
internal energy, the latter can be used as another desired output. The third output, the mass 
flow rate, is kept as a function of the states. Summarizing, the outputs are defined as 


V = 


P 


u 

C 3 ^/(P(p,u) - P s )pC ve 


(3.1) 


The corresponding desired quantities are precomputed from the desired pressure, temperature 
and mass flow rate profiles as follows: 
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1. With P s and the prescribed output temperature profile, find the enthalpy profile. 

2. With the enthalpy and mixer pressure profiles, find the internal energy and density pro- 
files. 

3. The desired mass flow rate profile is used directly. 

3.2 Relative Degree and Input Integration 

The model of Eqs. 2.14 along with the outputs defined above can be compactly written as 

P = fi(P(p,u),p)C vl + f 2 (P(p,u),p)C vg + f 3 (P(p,u),p)C ve (3.2) 

u = gi(P(p,u),p,u)C v i + g 2 (P(p,u),p, u)C vg + g 3 (P(p, u),p,u)C ve (3.3) 

y = [p u -Vf 3 {P{p,u),p)C ve ] T (3.4) 

where the function definitions are chosen to match Eqs. (2.14) and (3.1) and V is the constant 

mixer volume. The model has certain characteristic which prevents direct application of in- 
put/output linearization theory. It is seen that the third component of y is statically related 
to one of the control inputs, namely C ve . This implies that the partial relative degree of y 3 is 
zero. Conceivably, one could use a time-explicit control of the form 

C ve {t,p,u) — — - - 

Vf 3 (P(p, u),p ) 

to attain perfect tracking of the output mass flow rate and use the two remaining controls to 
force the flow to have the desired pressure and temperature. This approach suffers from two 
drawbacks: the most serious one lies in the fact that the resulting reduced two-input, two-output 
system may be rendered non-minimum phase by the choice of C ve (t, p,u). Even if the reduced 
system is minimum-phase, the lack of robustness of a feedforward law is still questionable. 

One way to get around this problem altogether is to augment the exit valve channel with 
an integrator, that is, let: 

Cye — U 

where v is a new control input. Now C ve is regarded as a state, and the resulting system is of 
third order, with three inputs and outputs. If the arguments of functions gi and /* are dropped 
from the notation, the new system equations become: 


flC v l + hCvg + fzPx >e 

(3.5) 

P 1 f 'r/ T P'l ( ' vg T g 3 C ve 

(3.6) 

V 

(3.7) 

[p u -Vf 3 C ve ] T 

(3.8) 
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Upon differentiating the outputs once, it is seen that the partial relative degrees are all 1: 


9 1 

92 

m 


flC v l + f'2C v g + hC ve 

9lCv i 92Cvg (]:*,( 'v<- 


( 3 . 9 ) 

( 3 . 10 ) 


= -V 


v f3 + C v 


P feCvg + f3C ve )+ 


d f OP 

+ ~Qp-Q^{9iCvi 92C vg + gsC ve ) 

+ d ^{hC vl +f 2 C v g + hC ve ) 

Upon rearranging, the output derivatives can be expressed compactly as 


( 3 - 11 ) 


’ ij\ " 


fsCve 

92 

= 

93 Pve 

_ 93 _ 


_ ~ VC leVw{%f3 + ^93) + ^f/3] . 



fl 


+ 

91 


h 

92 


VC ve [Jp ( |p/l + + g^fl] VC ve [Jp (|^/2 + |^52) + 3^/2] U/3 




dPl dp- 


dp 


or, 


( 3 . 12 ) 


y = D + Ew 

Provided E is invertible in a region U of the state space, the feedback law 

w = E~ 1 (y d - F(y - y d ) - D ) 
achieves exact linearization of the system, with tracking error dynamics given by 

(y - 9 d) + r(y - y d ) = 0 

If r is chosen as a diagonal positive-definite matrix, the resulting control law is called “decou- 
pling control”, since the dynamics of the outputs are decoupled. If the appropriate function 
definitions are substituted, the forms of the E and D matrices are as follows: 


fi = 
f2 = 

f2 = 

h = 

9i = 
93 = 


CW(Pi-P)pi 

C’ 2l \iP g >2P 

C'iy/p* - P 2 , if Pg < 2P 

cW(P-Ps)p 

cW(Pi-P)pi 



D = 


P-P S P 


p p 

cW(p - Ps)pCv 


n» 

°3 


VCleCj 


p-p a p r 

~ C ' ve 


P% + 0-185ff + (P-P S 


p du 
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cW( p i - p )pi 

E= gjVW - p )pi {*p) 

VC ve C{C^(P l -P)p l r gp , , t , BP , „ N 1 

* + ( h i -“)fe + ( p - P =)J - 2 vtp 


2 V(P~P S )P 


VC ve f2C £ 
Ps)p 


0 
0 

+ (fr’fl — u ) + (-P — P a )] ~V (P — P s )p 


where the constants C(, C 2 , C 3 and C 4 are related to those of Eq. 2.14 by C[ = C\/V, 
C' 2 = Ci^fTgPglV , C' 3 = -C 3 /V, C’i = 0.185 C’ 3 and C 4 = C 4y /T gPg /V. 


3.3 Input-Output Linearizability of Augmented Model 

The ability to construct a feedback linearization controller hinges, firstly, on stable zero dy- 
namics of the augmented system, and, secondly, on the invertibility of matrix E. 

3.3.1 Invertibility of E 

By inspection, it is readily seen that the first two rows of E are linearly independent provided 
h g 7 ^ hi. This has a direct physical interpretation: if the two fluids have the same thermal 
properties (i.e., enthalpies), the ability to change the thermal properties of the mixture by 
changing the relative flows is lost. Gas and liquid enthalpies are different for the expected 
operating conditions. The third row is linearly independent from the first two if P — P s > 0, 
which is also true for the mixer. Therefore E is invertible over the whole range of expected 
mixer operating conditions. 


3.3.2 Zero Dynamics of the Augmented Model 


As it is known, the zero dynamics is preserved under input transformations. This implies that 
one may examine the non-augmented model for zero dynamics and draw conclusions about the 
augmented model’s zero dynamics. Suppose the exit flow is to be held constant at a value Y 30 . 
The only way in which this can be achieved is by letting 

c “ (t) = ~vm (3 - 13) 

at all times. Differentiating the other two outputs and equating them to zero results in 

C v i(t) = (-f 2 ( t )C vg (t) + ^j (3.14) 

Cv9it) = -gtf) { 9l ^ Cvl ^ + 9 Vf!(t°) (3 ' 15) 

Upon substitution and rearrangement, it is seen that there exist three uniquely defined control 
inputs which hold the outputs constant. Since two of the outputs coincide with the states, it 
trivially follows that they are kept bounded, and therefore the system has stable zero dynamics 
(i.e., the system is minimum-phase). The control inputs are given by Eq. (3.13) and, dropping 
the time notation, 


c, 


v g 


C vi 


Do/l ( 93 _ 9l\ 

y y/3 h) 

f \ 92 - gif-2 

fjo . /1/2 ( 53 _ 9 i \ 

V h . fi92 - gif -2 \h fi) 


(3.16) 

(3.17) 
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Note that the above formulas can be used to find the valve positions at which the system has 
a prescribed outflow and a pair of thermodynamic properties. Therefore it produces the same 
results as the program mxr-cv described earlier. 


3.4 Implementation Issues 


Aside from the computational burden introduced by the inversion of E at each time step, 
the control law requires knowledge of the partial derivatives and ^ at all times. This 
significantly increases the complexity of the controller, since not only P(p , u ) is required as a 
table, but two extra tables are required which contain the partial derivatives. The complexity 
is somewhat reduced by noticing that one of such partial derivatives has a very small value 
for the range of properties of Hydrogen considered. In fact, all thermodynamic data in use for 
mixer simulations consists of a table of densities and internal energies as a function of pressure 
and temperature. That is, 

P = P(P,T ) 

and 

u = u(P, T) 

The data for Hydrogen between 2000 and 13500 psia and -400 and 200°T 1 reveals that 


du 

df 


» 


dp 

dT 


p 


(3.18) 


For the above pressure and temperature ranges, it is possible to show from the data that given 
density and internal energy, one can recover the pressure and temperature. That is, 


P = P{p,u) 


and 

T = T{p,u) 

are functions of two variables. Noting that P and T are independent variables we have 

dP _ dP dp dP 9u_ 
dT ~ ~dpdT + ~du&T ~ 1 

and making use of Eq. 3.18 it is seen that 

dP 

~df~° 

Thus the above partial derivative is set to zero in matrices D and E. The other partial derivative 
can be found as 

dP _ 1 

The derivative on the right is readily found from the data matrices. 

Even though it is unlikely that the control resources available can accommodate even one 
partial derivative, a simulation study using derivative information is carried out next. It pro- 
vides a basis for comparison of other designs, since the feedback linearization controller with 
full information is expected to give excellent performance. 


30 


RELEASED - Printed documents may be obsolete; validate prior to use. 



Parameter 

Value 

Units 

Pi 

8500 

psia 

Ti 

-340 

°F 

P* 

13500 

psia 

T a 

90 

°F 

Ps 

5533 

psia 

Po 

3.85 

C? 

CO 

U 0 

180 

Btu/lbm 


Table 3.1: Feedback Linearization Parameters 

3.5 Simulation Study of the Feedback Linearization Controller 

The purpose of the study is to test the suitability of the partial derivative approximations and 
to provide insight about the limits to attainable performance. Initial simulations will hold plant 
parameters constant and assume that valves open instantaneously. Then the thermodynamic 
properties of the supply fluids will be perturbed to reflect actual operating conditions. First 
order dynamics and saturation will be imposed on the valves and the performance of the nominal 
controller evaluated. 

3.5.1 Nominal performance 

Table 3.1 summarizes the parameters used in the simulations. The decoupled control gains were 
set as 10 for the density, 10 for the energy and 5 for the exit flow. The controller has exact 
knowledge of these parameters, except the initial conditions. The objective in Simulation 1 is 
to transfer the mixer from the initial conditions to a pressure of 6000 psia, exit temperature of 
— 200°-F and exit flow of 40 lbrn/s. With the known and fixed back pressure it is determined 
that the mixer states must be transferred to p = 2.93 Ibm/ ft 3 and u = 400.9 Btu/lbm for the 
pressure and exit temperature to converge to their desired values. The reference states and 
flow are specified as constants, leaving the trajectories to be determined by the controller. In 
Simulation 2, the mixer is asked to reach the previous equilibrium point and then raise the exit 
temperature to 0 °F and the pressure to 7000 psia, while reducing the exit flow to 10 lbm/s. 
This corresponds to lowering the density to p = 2.143 and raising the energy to u = 995.9. This 
time, the states and the flow are required to track ramp functions which transfer from old to 
new values in 5 seconds. Figs. 3.1 and 3.2 show the results of Simulations 1 and 2, respectively. 
The steady values required are exactly achieved. However, the transient response of the exit 
flow in the simulations has overshoot. This may be reduced by modifying the control gains, 
but perhaps at the expense of slowing down the response. 
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Figure 3.1: Simulation of Feedback Linearization Controller 
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Figure 3.2: Simulation of Feedback Linearization Controller 
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Chapter 4 


Sliding Mode Control of the Mixer 
with Valve Dynamics 

4.1 Mixer Model with Valve Dynamics 

The developments of the previous sections assume that there is direct and instantaneous com- 
mand over the valve opening coefficients, which are taken as control inputs. In reality, a C v 
coefficient is the response of the valve to a given command. The valves have their own analog 
controller, which makes the valve behave as a first order linear system with certain gain and 
time constant. Thus, it is assumed that the valve models are given by 

C v = A c C v + B c C (4.1) 

where A c = diag(— 1 / ti , — 1 / 72 , — 1 / 73 ), B c = [fa fa fa] T ■, T i are the time constants of the 
liquid, gas and exit valves and fa are their respective gains, Q = [£i Q 2 Cs] T is the vector of 
control inputs to the liquid, gas and exit valves, and C v = [C v i C vg C ve ] T is the vector of valve 
opening coefficients. The C\, coefficients now become states in the model equations, which can 
be written as 


p = 

1 T 
v c 'f 

(4.2) 

u = 

\c T v A p f 

(4.3) 

Cv = 

A c C v + B c ( 

(4.4) 

V = 

[p U C ve fe] T 

(4.5) 


where / = [f\ /2 — f e ] T and A p = diag(^p^, ha p u , 0.185^). The feedback linearization 
controller presented earlier has several disadvantages which limit its applicability in the real 
system. One of the most salient disadvantages is its lack of robustness due to dependence on 
an accurate model and constancy of parameters. Both can be handled by means of Sliding 
Mode Control, a technique which is specifically tailored to work with multivariable systems 
with certain types of disturbances and parameter variations. In the specific case of the mixer, 
the most significant unmodeled effect is the heat transfer with the environment. Although 
possible to model in principle, a quantitative description of the heat transfer through mixer 
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walls involves awkward convection correlations with parameters which are heavily dependent 
on temperature and flow characteristics. The convection coefficients can be extremely uncertain 
and often have to be found through iteration, even for the steady heat transfer calculations. 
For control purposes, it is far more convenient to treat the heat transfer (and other unmodeled 
energy losses) as a lumped and unknown bounded disturbance. A bound on the external 
heat transfer is relatively easy to find. As it will be shown in the next sections, the heat 
transfer appears in the equations as a matched disturbance, to which the system can be made 
insensitive through Sliding Mode Control. The feedback linearization controller of the previous 
sections was derived with the implicit assumption that the pressure and temperature of the 
fluids reaching the mixer are constant. While a feedback linearization controller can be derived 
that takes into account those variations, it would require knowledge of the inlet thermodynamic 
properties at all times, and what is worse, the derivatives of the properties would also be 
required. Performance tradeoffs and limitations are expected, however, due to the intrinsic 
limitations of valve actuators in the system. 


4.2 A Pressure-Density Model with Heat Transfer Disturbance: 
Ideal Gases 


It is convenient to start the study of Sliding Mode controllers by simplifying the mixer model. 
The model will be simplified by assuming that the fluid reaching the inlet valves is an ideal gas, 
which obeys the well-known state equation 


P = pRT 


where R is the ideal gas constant in ft 3 -lbf/in 2 -lbm-R, T is the temperature in degrees Rankine, 
P is the pressure in psi and p is the density in lbrn/ft 3 , for English units. For SI units, the 
choices are also straightforward. It is convenient to write the equations in term of a pressure 
state instead of internal energy, since the flow rates can be directly computed from the state 
variables without further reference to thermodynamic data. For the case of the ideal gas, the 
following identities and definitions are in order 


h = 

C T - PCp 
p pR 

(4.6) 

u = 

CT - PCv 
v pR 

(4.7) 

k = 

Cp 

C~V 

(4.8) 

R = 

Cp-C v 

(4.9) 


where C p and C v are the constant pressure and constant volume heat capacities, respectively, 
and k is the heat capacity ratio. Using these definitions and identities, and performing manip- 
ulations, the A p matrix of Eq. (4.3) can be written as 


A 


p 


diag 


Rh\ 

gT’ 


Rh>2 kP 

~gT’ T 


(4.10) 
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The A p matrix is dependent on the state, as it can be seen in the above formula. The heat 
transfer disturbance is directly included in the basic energy balance of Eq. 2.4, where the 
subindices have been modified to reflect the ideal gas assumption: 


mu + um = m\h\ + 777 , 2/12 — fh e /i e + Q 


where Q denotes the heat transfer rate, expressed in Btu/sec or in Watts, according to the 
system of units. The energy balance equation can be written in terms of the pressure derivative 
as 


P 


1 

V 


Rh\ Rhv P 

— 7771 + -^—1712 - k—m e + 
(-Sy Csy P 


RQ\ 

C v ) 


Using the above definitions and identities, the final form of the ideal gas mixer model with heat 
transfer disturbance is given by 


p = 

v c ' f 

(4.11) 

ii = 

l C V A pf+ y c Q 

(4.12) 

Cv = 

A c C v + B C Q 

(4.13) 

y = 

[p U C ve f e ] T 

(4.14) 


with A p given by Eq. 4.10. 


4.2.1 Equilibrium Points 

It is convenient to start the study of model properties by considering its equilibrium points. It 
is of particular interest to characterize the equilibrium points which correspond to positive flows 
and positive valve openings. This information is of crucial importance in the operation of the 
mixer, for the model may allow for operating conditions which are not of practical interest, such 
as those necessitating reverse flows. Note that the equilibrium analysis alone is not sufficient 
to guarantee that a particular operating condition is reachable with given system resources and 
without violating physical constraints during transient conditions. 


4.2.2 Sliding Mode Control Design 


Having established the stability of the zero dynamics, it is now possible to follow the steps of 
a standard sliding mode controller design. The first step is to write the 5th-order system in 
terms of output derivatives. Differentiating the first output twice leads to 


Cl 


V 


yi - ~rr + -77 + 


df\ , C T B c f 


dt 


V 


where 


dj_ 

dt 


dfi dh 
dt dt 


dfe 

dt 


Similarly, differentiating the second output twice gives 

dp ' T 


V2 = 


V 


a a dAp 

AcAp H P~ 

dt 


, + Ap Tt 


C ‘_BcArj 

V 


R 

VC, 


Q 
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The third output, however, needs to be differentiating only once for the control to appear: 


m = c ve d ^ + (--c ve + p 3 ( 3 )f e 

at 73 

The sliding manifolds are chosen to yield linear tracking errors when the sliding regime has 
been reached. Specifically, choose 


•Si 

= Pd- pP c p (p d - p) 

(4.15) 

S2 

= P d -P + cp(P d -P) 

(4.16) 

S3 

= C w (wd C ve fe) 

(4.17) 


where p d , Pd and Wd denote the desired trajectories for the density, pressure and exit mass flow 
outputs, respectively, and c p , cp and c w are sliding coefficients. It is desired that the si, S 2 and 
S 3 reach zero after some finite time and remain zero despite of the heat transfer disturbance. 
Choosing Lyapunov functions 



for z = 1,2,3, it is sufficient to ensure that V, < 0 at all times to guarantee the convergence of 
Si . To reach zero in finite time, a commonly used choice is to enforce 


V = SiSi < - rn\si\ 


(4.18) 


It can be shown [30] that the above implies that s, reaches zero in a time ti given by 

\si(t = 0 )| 


tj < 


Vi 


Now, the derivatives of the sliding functions can be obtained as 


• _ , /• 1 nT f \ C T Bcf C% ( df 

Si - Pd + c p \ Pd - yC v f J — — lA c f + — 


• p 

S 2 ~ Pd-Py 


A A + \ f + a — 

AcAp+ dt r +Ap dt 


dt 

C T B c A p f 

V 


R 

VC, 


Q + cp ( Pd — —C^Apf — 


RQ 

VC v 


£3 — C-u 


Wd - ( ~^C ve + /3 3 Cs) fe ~ 


Enforcement of Eq. (4.18) can be achieved by using the equality for si and S 3 , since the heat 
transfer uncertainty is not involved in the expressions, that is, let SjS, = —r]i\si\ for i = 1 , 2 . 
This is equivalent to Si = — r/isign(sj). For S 2 the equality cannot be guaranteed, since it 
would require knowledge of Q and Q. The Sliding Mode Controller will rely, however, on the 
knowledge of a bound for these quantities and on the ability to increase the sliding gain r / 2 to 
counteract the disturbance. Enforcing s* = — r)isign(si) for z = 1,3 leads to 


C T Bcf 

V 

C T cf 


id 

r 3 


(4.19) 

(4.20) 
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where 


C T ( df\ 

hi = Pd + CpPd + r/isign(si) - -y- f (A c + c p I)f + — J 

r 3 = c w + ~^C ve f e - Cve^f) + %sign(s 3 ) 

and G is a 3-by-3 matrix with zeros in every entry, except for G( 3, 3) = — c^/3 3 . For the 
remaining sliding variable, we choose 


where 


so that 


( T B c A p f 

V 


C T 

^2 = Pd + c P P d + 7/ 2 sign(s 2 ) - -y- 


= r 2 


(4.21) 


dA r 


df 


S 2 = -?? 2 sign(s 2 ) - 


(A c + c P I)A p + ^R\f + A p J 


RQ c p R 


We want 
with rj 2 > 0. Let 

where Q is a known bound such that 


VC v VC v 
S2S2 < -V2\ s 2 


, ^ RQ 

m = v 2 + 


(4.22) 


VC' 


\Q + cpQ| < Q 


at all times. Then it is straightforward to show that Eq. 
thus sufficient to choose 


m > 


RQ 

VC v 


(4.22) is satisfied. Since rf 2 > 0, it is 

(4.23) 


. The equations that yield the control input, Eqs. (4.19), (4.21) and (4.20) can be compactly 
expressed in matrix form as a system of equations which is linear in Q. 


1 

0 

1 


' Fi ■ 

f T A p B c 

V 

c = 

r 2 

f T G 


. r 3 . 


(4.24) 


Invertibility of the system matrix is related to the zero dynamics and to the existence of feasible 
trajectories, which is currently under study. 
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Parameter 

Value 

Units 

Pi 

6.894 

MPa 

Ti 

311.1 

°K 

P2 

5.516 

MPa 

t 2 

422.22 

°I< 

Ps 

206.842 

kPa 

R 

296.8 

J/kg°iF 

c p 

1041.6 

J/kg°AT 

k 

1.398 

n.a 

Po 

3.447 

MPa 

Po 

48.05 

kg/rn 3 

C v 1,2,3 

10 

n.a 

71,2,3 

0.1 

sec 

Pi 

1 

n.a 

th 

2 

n.a 

P3 

3 

n.a 

V 

0.07079 

m 

c 

2.40428 xlO -5 

n.a 


Table 4.1: System Parameters for Sliding Mode Controller Simulation 

4.2.3 Simulation Example 

In this example, we choose the SI system of units with Nitrogen as the working gas and the 
following design specifications: 

• Obtain a steady exit flow of 10 kg/s with a mixer operating pressure of 4 MPa. 

• The temperature of the exit flow must be 337 °K = 64 °C. 

• The heat transfer disturbance is is Q = 10sin(10t). 

• The steady operating point must be reached in less than 25 seconds, for known initial 
conditions given in Table 4.1. 

• Valve chattering must not be observed. 

Other system parameters are summarized in Table 4.1. The value of c in Table 4.1 is the 
appropriate one for the valve flow formula in SI units. The formula for liquids has been chosen 
for simplicity and illustration purposes. 

Sliding Mode Controller Design The values of w ( i and P,i are directly taken from the 
design specifications. It is necessary to calculate the value of pd which will result in the desired 
exit temperature. For this, we have that the isoenthalpic process at the exit valve is also 
isothermal, since it is assumed that C p is constant. Therefore we choose pd = = 40kg/m 3 . 

In order to meet the settling time requirement, we choose c p = cp = c w = 100, rji = rjs = 120 
and i ]2 = 2.4062 x 10 6 . With the known initial conditions and setting Pd = Pd = Pd = Pd = 
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Wd = 0 we obtain the reaching times: 


tr 1 

< 

tr2 

< 

tr3 

< 


\si{t = 0 )| 

m 

= 0)1 
m 

M* = o)| 

% 


7. 08-sec 
19. 5-sec 
5. 83-sec 


The overall settling time will be given by the sliding surface reaching time plus the settling time 
during sliding, which is given by t s \ = t s 2 = t s 3 = 4/100 = 0.04sec.. The above gain choices 
thus satisfy the settling time requirements. Finally, we check that the controller can attain its 
objective despite the presence of the heat transfer disturbance. The value of Q is computed as 


Q = \Q + c P Q\ = 1100 


We see that the selection of ij 2 is appropriate, since 

QR 

i n = 2.4062 x 10 6 > 

V U v 

To avoid chattering, the signurn functions present in the control law are replaced by finite slope 
approximations consisting of saturation functions: 

sign(sj) « sat (si/fa) 


where 


sat(x) 


x if |a:| < 1 
< 1 if x > 1 
—1 otherwise 


(4.25) 


For the simulation example, the values (f>i = <j >3 = 10 and (f > 2 = 1 x 10 4 were chosen. Fig- 
ures 4.1, 4.2 and 4.3 show the simulation results. As it can be seen, the control objectives are 
satisfied. The second valve position history, however, is physically meaningless, since it contains 
a time interval when the valve opening is negative. Avoidance of reverse flows and negative 
valve positions was not specified in the control problem. This simulation motivates for further 
study into admissible trajectories and controllers which produce them. 
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Valve 1 Coeff. Valve 2 Coeff. 




Exit Valve Coeff. Heat T ransfer Disturbance 




Figure 4.3: Simulation of Sliding Mode Controller - Control Variables and Disturbance 
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Chapter 5 


Constrained Sliding Mode Control 
Design Using Robustly Positively 
Invariant Cylinders 

5.1 Introduction 

A practical solution to the robust control problem of the mixer under constraints requires an 
examination of invariant set theory. Valve positions are physically restricted to be nonnegative 
and bounded by certain constant value related to allowable stem travel. In the model derived 
above, such positions are state variables, while the commanded positions are control inputs. 
It is therefore important to examine conditions and design methods which guarantee that the 
controlled states belong to certain predetermined operating set, even under uncertainty and 
disturbances. In mathematical control theory, such operating sets are termed invariant, and 
several kinds of invariance are defined. A comprehensive survey of set invariance in control 
theory is out of the scope of this report, however some definitions will be required. For details, 
please refer to [4, 8, 20]. Much of the work in set invariance concerns linear systems and 
ellipsoidal or polyhedral operating sets. New theory is developed and introduced here for 
linear and a class of nonlinear systems under Sliding Mode Control. The invariant sets under 
consideration here are hyper-cylinders which project orthogonally as ellipsoids in a reduced 
dimension space. 

5.2 Basic Theory of Robust Positive Invariance 

Although the systems to which the main results of this article apply are linear, the results of 
this section apply to more general systems, possibly time-varying. The explicit dependence 
on time has been dropped from the notation, however. Given a dynamic system described by 
x = f(x(t)), having a unique solution x(t) in a subset fl C M n , the set S C is said to be 
positively invariant (PI) for the system if for every initial state x Q £ S , the solution x(t) belongs 
to S for t > 0. For the uncertain system x = f(x(t),w(t)) where w(t) is an input known to 
belong to some set W, S is said to be robustly positively invariant (RPI) for the system if for all 
x 0 £ S and for all w(t) £ W the solution x(t) belongs to S for t > 0. The widely known result 
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due to Nagumo [17] provides a necessary and sufficient condition for the invariance of a set S 
in terms of its tangent cone, also known as Bouligand contingent cone. A precise definition of 
tangent cone is given, for instance, in [4] and the reader is referred to it for details. For our 
purposes, it will suffice to state what the tangent cones of particular sets are, and to cite some 
relevant properties. The notation Cs(x) will be used for the tangent cone of set S at a point 
x £ R n . The notations A c , A, dA and A C B indicate, as usual, the complement, closure and 
boundary of A, and that A is a subset (not necessarily proper) of B, respectively. The extended 
real line is denoted here as R* = M U { 00 } U {— 00 }. 

Theorem 1 (Nagumo [17]). Consider the system x = f(x(t)) having a unique solution for 
each initial condition in a set Let S C 11 be a closed and convex set. Then S is PI for the 
system if and only if f(x) £ Cs(x) for all x £ dS. 

Mechanically interpreted, the theorem formalizes the intuitive notion of invariance being 
attained when the velocity vectors x at the boundary of the set all point into, or are tangent 
to S. The following important extension of Nagumo’s result is given in [1] and concerns robust 
positive invariance. 

Theorem 2. Consider the uncertain system x = f(x(t),w(t)) where the uncertain input w(t) 
has values in W for all t > 0. Assume that the system posesses a unique solution for all initial 
conditions xq £ 11 and all w(t) £ W. Then the convex and closed set S C Ll is RPI if and only 
if f(x, w) £ Cs(x) for all x £ dS and for all w £ W. 

A few useful properties of tangent cones are now stated. 

Property 5.2.1 ( [1, 12]). If A and B are closed and convex sets such that 0 £ int(A — B), 
then Cac\b{x) = Ca(x ) n Cb(x ) for all x £ An B. 

Property 5.2.2 ( [1, 4]). If A is convex, then Ca(x ) is convex for all x £ A. 

It is also a fact [1, 4] that if x £ int(A) then Ca(x) = M n . To end this section, we show 
that, for certain class of systems, and when the uncertain input belongs to a closed interval, it 
is sufficient to check for positive invariance of the system at extreme values of the input. 

Theorem 3. Consider a system that is linear in the uncertainty, e.g., x = f(x(t),w(t)) = 
g(x(t)) + Bw(t) where w(t) £ T = [w , w\. Assume that the system posesses a unique solution 
for all initial conditions xo € fl and all w(t) £ T. Then S C fi is RPI if and only if it is PI for 
x = f(x(t ) , w) and for x = f(x(t ) , w) . 

Proof. Suppose S is positively invariant for x = f(x(t),w) and for f(x(t),w). Then, by Theo- 
rem 1, f(x,w) £ Cs(x) and f(x,w ) £ Cs(x) for all x £ 8S. Let w £ 1. Then w = aw + (1 — a)w 
for some a £ [0, 1] and thus f(x, w) = af(x, w) + (l — a)f(x, w). Since Cs(x) is convex whenever 
S is convex, we conclude that f(x, w) £ Cs(x) for all x £ 8S and for all w £ 2. By Theorem 2, 
S is RPI. The reverse implication is trivial. □ 


5.3 Sliding Mode Control of Linear Systems 

Sliding Mode Control (SMC) is a widely studied technique [10, 30, 31] that achieves total insen- 
sitivity of the controlled variables to certain kinds of disturbances and parameter uncertainties. 
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In this section we briefly introduce the salient characteristics of linear systems under SMC. 
Consider the single-input, linear and time- invariant system 

x = Ax + Bu(t ) + D5(t) (5-1) 

where A is an n-by-n matrix, B and D are column vectors, x 6 R n , u is the scalar control input 
and 5{t) is a scalar, unknown disturbance. It is assumed that the pair (A, B ) is controllable 
and that the matching condition 


rank[S|D] = rank(L>) 

is satisfied, so that the system can be rewritten as 

x = Ax + B(u(t) + C(t)) (5-2) 

Assume that S(t) is bounded so that ((t) e Z for all t > 0, where Z = [— (, £]. In order to 
specify a sliding mode control law with linear sliding manifold, let T be a nonsingular matrix 
such that T~ 1 B = [0 0...| 1] T and such that An is a nonsingular n — 1 by n — 1 matrix, where Ajj 
for i,j = 1,2 are the partition blocks of T~ 1 AT such that A 22 is scalar. Such transformation 
can always be found. Note that taking T as the transformation matrix for the control canonical 
form is not a good choice, since in that case the first column of An is zero. A method to find 
such a transformation is given in the appendix. Define a coordinate transformation x = Tz 
and write the system equations in the new coordinates as 

z = T~ l ATz + [0 | 1 } T (u{t) + C (*)) 

Consider the sliding manifold s = Gx = GTz = [G'i | G^z. Without loss of generality, consider 
that G 2 = 1. The control law 

u(t) = -A 21 Z 1 - A 22 Z 2 - Gi(An 2 :i + A 12 z 2 ) - 11 sgn(s) (5.3) 

results in the closed-loop dynamics described by 

zi = A u zi + A U Z 2 

z 2 = -GiAnZi - G 1 A 12 Z 2 - r 1 sign(s) + C(t) 

S = G\Z\ + Z2 

It can be easily shown any choice of G under the above constraints on T, and such that G 2 = 1, 
results in GB = 1. It is likewise straightforward to show that the closed- loop dynamics in the 
original coordinates is described by 

x = (I — BG)Ax — Br) sign(s) + B((t) (5.4) 

Using V = s 2 as a Lyapunov function shows that if i] > ( then control law (5.3) results in 
the state reaching the plane s = 0 in finite time and remaining there indefinitely despite the 
presence of the disturbance [10, 30, 31]. Evolution of the closed-loop system (5.4) for s = 0 is 
independent of the disturbance and described by the reduced dynamics 

zi = A w z 1 (5.5) 
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where 


A w — Ai ~ A 2 G 1 (5.6) 

Thus, a stable sliding mode is obtained by choosing G\ such that A w has eigenvalues with 
negative real parts (A w is Hurwitz). It can be shown that the controllability of ( A , B) guarantees 
that the eigenvalues of A w may be freely placed using G\ . For the remainder of this article, it 
will be assumed that A w is Hurwitz, 7 > ( and that a unique solution to the closed-loop SMC 
differential equation (5.4) exists for every initial condition in R n . 

5.3.1 A Useful Decomposition of the Closed-Loop Dynamics 

The Lyapunov function V = s 2 induces an obvious family of invariant sets, namely the sets 
S 7 = {z € M n : |s(z)| < 7} for 7 > 0 are all positively invariant. These sets are “naturally” 
invariant for systems under SMC. A coordinate transformation is introduced here that decouples 
the motion towards s = 0 from the overall convergence to the origin. This decomposition will 
suggest a cylindrical shape for positively invariant sets. 

Lemma 5.3.1. There exists a coordinate transformation x = Jw with J nonsingular in which 
the closed-loop dynamics (5-4) is expressed as 



where B w = [Bi\B 2 ] 7 = — J 1 B and G w = — G\A^ A 12 + 1. Moreover, B 2 G W = —1. 

The proof to Lemma 5.3.1 is done directly by specifying J and is shown in the appendix. 

Note that s is just proportional to the scalar 1 V 2 and independent of w \ . An immediate obser- 
vation is that an arbitrary (possibly infinite) real interval is RPI for the dynamics of w 2 . The 
result is formalized in the following 

Lemma 5.3.2. For any initial condition [rcf 0 | W 2 o] T £ R n , where W 20 belongs to the possibly 
infinite interval [a,b\ 3 0, the trajectory w(t) of the closed-loop system (5.7) is such that for all 
W 2 (i) £ [a, b] for all t > 0, that is, [a, b] is RPI for the dynamics of W 2 - 

Lemma 5.3.2 follows directly from Lemma 5.3.1 by using the conditions B 2 (—G\Af^A \2 + 

1 ) = — 1 < 0 and 7 > Q in establishing the monotonic decrease of |s| (and that of |w^|) at 
both sides of s = 0. With the aid of Lemma 5.3.2 it is now possible to specify cylindrical 
positive-invariant sets with fairly general cross sections. Introduce the notation 

TLs 0 [a,b\ = {[wJ\w 2 ] T £ K n : wi £ S Q , w 2 £ [a,b]} 

The set S Q is termed the cylinder’s cross section. The following result follows directly from 
Lemma 5.3.2: 

Theorem 4. Let S a C M n_1 be a compact and convex set containing the origin . Suppose S Q is 
RPI for the system w \ = A w w\ + Bi<4'(t), where A w is Hurwitz and |£ 7 (£) | < 7 + ( for all t > 0. 

Then all cylinders TLs 0 [a, b] such that 0 £ [a, b] are RPI for the closed-loop dynamics of (5.7). 

From now on, a compact and convex cross section will be assumed in the notation for cylin- 
ders. Also, the notation Hs„ is used for the family of cylinders TLs 0 = {Hs 0 [a,b] : 0 £ [a, b], a,b £ M*}. 
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5.4 State- Constrained Invariant Sets 


Robust positive invariance of a set S alone is not sufficient to guarantee that the state of the 
controlled system will not exceed prescribed bounds given by a set G (except in the uninteresting 
case S C G). Consideration here is restricted to linear state constraints. 

Definition 5.4.1. A linear state constraint set is defined as G = fl Gi for i = 1,2. .to, where 
Gi = {w £ M n : Tju; < 1}. T* are row vectors such that G is a convex set containing the origin 
in its interior. 

Note that 8Gi , the boundary of Gi, is the set of points such that rite = 1. In the context 
of SMC under state constraints, G is assumed to be the constraint set in w coordinates. That 
is, to each linear constraint Ti x x < 1 for the original system there corresponds a constraint 
Ti x Jx = TiW < 1. Since J is nonsingular, it is straightforward to see that G is a convex if and 
only if the original constraint set in x coordinates is so. Furthermore, assume that T i x B 0 
for i = 1,2.. to. This basic feasibility assumption for the problem goes beyond controllability of 
the pair {A, B) [9]. At this point it is convenient to state that the tangent cone to a linear state 
constraint Tjtc < 1 is given by TiW < 0. In connection with the uncertain system w = f(w, £(i)), 
where £ £ Z for all t > 0, introduce the following sets: 

Gf = {re € M" : Tif(w, () > 0 for some ( £ Z} 

G~ = {w £ M n : Tif(w, C) < 0 for all ( £ Z} 

(5.8) 

The next result concerns the invariance of the intersection of a state constraint set and an 
RPI set. It generalizes Theorem 1 from [20], in that the system is not restricted to be linear 
under state feedback, the positively invariant set does not have to be an ellipsoid, disturbance 
is allowed and the constraints may be asymmetric. 

Theorem 5. Let M C R n be a compact and convex set which is RPI for a system w = 
f(w, where £(i) £ Z for t > 0. Let G = n Gi, i = 1, 2. .to be a linear state constraint set as 

in Def. 5-4.1. Assume, further, that 0 £ int(M — G). Denote S = M n G. Then S is robustly 
positively invariant for the system if and only if S n 8Gi fl Gif = 0 for i = 1, 2 ...to. 

The proof of Th. 5 is shown in the appendix. Since G is compact, the following quantities 
are well-defined: 


w 2 = max W ‘2 (5.9) 

w£G 

w 2 = m in W 2 (5.10) 

w£G 

Note that w 2 < 0 < W 2 , since 0 eint(G). The following Corollary provides a sufficient condition 
for robust positive invariance that leads to computation. 

Corollary 5.4.1. Let the family of cylinders 7is 0 be RPI for the SMC dynamics (5.7). Let 
G be a constraint set such that 0 £int(G — TLsf)- Then S = TLs 0 H G is RPI for (5.7) if 
Hs 0 [ w 2 ,W2\ n dGi fl dGf = 0 Vi = 1,2 . ..to. 
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For the remainder of the development consider 

f(w, 0 = [A w w i + B 1 (r]sign(G w w 2 ) - () \ B 2 (r]sign(G w w 2 ) - C)] T (5.11) 

Eq. (5.11) defines as the value of w in the SMC closed-loop dynamics for a constant 

disturbance value. The following technical results are useful in the proof of Corollary 5.4.1: 

Lemma 5.4.1. Let the family of cylinders 7is a be RPI for the SMC dynamics (5.7). Let a 
linear state constraint Gi be such that dGi n Bs 0 [w 2 ,w 2 ] / 0. Then Eire £ Hs a [ ul 2 , w 2 \ D dGi 
such that Tif(w,C) < 0 VC £ Z 

Proposition 5.4.1. Suppose Tj = [Fji | Ti 2 ], where / 0. Then the collection of subsets 

{int(Gf),int(Gf),dGf} corresponding to f(w,() of Eq. (5.11) is a partition ofW 1 . Further- 
more, 

dCf = {weR n : = 0 ,w 2 + 0}u{ [wj | 0] T £ M n : |r a A w u;i - r t B w (*\ < \T t B w | ?? } 

where (* = — sign(F iB w )( . 

The proof of Lemma 5.4.1 and Proposition 5.4.1 are shown in the appendix. 

Proof of Corollary 5.f.l. By hypothesis, and using Proposition 5.4.1, we have that either 7is 0 [w_ 2 ,W2\C 
dGi Cint(G^) or Hs 0 [w 2 ,w 2 ] n dGi Cint(G“). For the first possibility we have that 

hLs 0 [w 2 ,w 2 ] n dGi C int(G+) C Gf 

Then, for all w £ Bs 0 [w 2 ,W 2 \CdGi , 3C £ Z such that Tif(w, () > 0, contradicting Lemma 5.4.1. 
Therefore it must be that TLs 0 [w 2 ,w 2 ] n dGi Cint(G“). But int(G“) n Gf = 0, therefore 
hi So [m 2 , Lo 2 ] fl dGi l~l Gf = 0. Finally, since S C hls 0 {— oo,oo) we have S n dGi H Gf = 0. The 
Corollary now follows from Th. 5. □ 

5.5 Cylinders with Ellipsoidal Cross Sections 

The results presented above apply to a fairly large class of cylinder cross sections S Q . Ellipsoidal 
invariant sets are extensively used in computations due to the correspondence to quadratic 
Lyapunov functions and their simplicity in relation to the also commonly used polyhedral sets. 
Ellipsoidal sets, however, can be conservative. The semi-ellipsoidal sets introduced by O’Dell 
for linear systems under state feedback achieve a good compromise between simplicity and 
conservativeness. A semi-ellipsoidal set is obtained by intersecting a linear state constraint set 
and an ellipsoidal invariant set [19, 20]. In this section we restrict S Q to be ellipsoidal and 
develop results leading to calculations for constrained SMC design. As seen in the previous 
sections, a cylinder is required to be RPI for it to yield an RPI set upon intersection with the 
state constraints. By Th. 4, we see that the cylinder cross section is required to be itself RPI 
for a linear system driven by a bounded disturbance. Moreover, by Th. 3, it is sufficient to 
establish invariance for extreme disturbance values. 
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5.5.1 RPI Cylinders 

Consider the system 

wi = A w w i + (5-12) 

where £ [— £']. We wish to find conditions under which the ellipsoidal set £ = {uq £ 
M n_1 : wj Pw\ < 1} is RPI for the dynamics (5.12). To this effect, note that [4] the tangent 
cone of £ at its boundary d£ is given by Cs(w i) = {y : 2 wj Py < 0}. Nagunro’s condition 
results in 

2wi P(A w w\ + Bi('(t)) < 0 along wj Pw\ = 1 (5.13) 

Conditions on P for (5.13) to hold can be derived following two approaches. One involves linear 
matrix inequalities (LMI) and it lends itself to ellipsoid volume optimization. The other ap- 
proach has a simpler form and provides the maximum sum of disturbance bound and switching 
gain alowable for a particular ellipsoid. 


5.5.2 LMI Approach 


Nagunro’s condition (5.13) is equivalent to a quadratic boundedness requirement, which, as 
shown in [5], can be equivalently expressed by the LMI 


P Ay, -f- A^jP -f- aP P B\ 

Q'BjP -al 


< 0 


(5.14) 


where P = P' > 0 is sought. 

Lemma 5.5.1. The set £ = {uq £ M n : wj Pw\ < 1} is RPI for system (5.12) there exist a 
symmetric, positive- definite matrix P and a scalar a > 0 such that the LMI (5.14) holds. 


Note that LMI (5.14) is always feasible, since A w is Hurwitz. Moreover, only values of a 
less than the maximum decay rate need to be considered, that is, a satisfies 0 < a < 2 |Re( A) | , 
where A is that eigenvalue of A w with the largest real part in absolute value [5]. The above 
matrix inequality can be readily solved using, for instance, the Matlab LMI Toolbox. 


5.5.3 Critical Switching Gain 

A computationally simpler alternative is obtained by using Theorem 3. Positive invariance is 
now required for the two autonomous systems that result when either —<fi or Q' is used. In view 
of the symmetry of P, Nagumo’s condition can be restated as 

wj (Ay,P + PAyfiwi ± 2('BfPwi < 0 (5.15) 

along wf Pw± = 1 


Define 

AlP + PA w = -Q 


It is straightforward to prove that a necessary condition for (5.15) to hold for some P > 
0 is that Q > 0. Thus we may assume that P is the unique symmetric, positive-definite 
solution to the above Lyapunov equation for arbitrary Q > 0, guaranteed to exist due to A w 
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being Hurwitz. Given Q = Q 1 > 0 and P = P' > 0, there exists some (3 > 0 such that 
Wi(A^P + PA w )w\ = —WiQw\ < — / 3wJ Pwi for all uq G d£. The quantity (3 is bounded by 
(3 =min wjQwi : uq G d£. The bound 0 is readily obtained by Lagrange multipliers and 
given by 



(5.16) 


where y 2 = v T Pv and v is that eigenvector of P~ l Q yielding the least value for the right 
hand side of Eq. (5.16). Knowing [3 it now follows that inequality (5.15) holds in d£ if 
±2£ / Bj P w\ < (3. Moreover, since (3 > 0, we consider 


2('\BjPwi\ < (3 

along tcf Pwi = 1 


(5.17) 


To obtain a condition equivalent to (5.17) we solve 

max 2C,’\Bi Pw\\ s.t. 

T' 

w l Pw\ = 1 


This is readily solved using Lagrange multipliers. The maximum is 2(, r Bf PB\, which must 

be less than (3. The resulting invariance condition can be also shown to be necessary, and it is 
summarized in the following 

Lemma 5.5.2. The set £ = {uq G M n : w~[ Pw\ < 1} is RPI for system (5.12) if and only if 
there exist symmetric, positive- definite matrices P and Q such that 


AlP + PA 


C'< 



— Q and 


with (3 defined in (5.16). 

Lemmas 5.5.2 and 5.5.1 lead to computation of invariant cylinders for SMC, as summarized 
in the next result. 


Theorem 6. Let £ = {uq G M n : wfPwi < 1}. The family of cylinders Tig is RPI for the 
SMC closed-loop dynamics of Eq. (5.7) if and only if there exist symmetric, positive- definite 
matrices Q and P and a scalar a > 0 such that Lemma 5.5.2 or Lemma 5.5.1 are satisfied, 
where (3 is defined in (5.16) and f = i] + (, with r) > (. 

Th. 6 is a direct consequence of Theorems 4 and 3. 


5.5.4 State-Constrained RPI Cylinders 

In this section we apply Corollary 5.4.1 to cylinders with ellipsoidal cross sections and provide 
results leading to design calculations. The shape of the set OGf and the partition it induces is 
schematically depicted in Fig. 5.1 for n = 3, as viewed from a plane perpendicular to the planes 
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Figure 5.1: The set dGf and the partition it induces. 


TnA w =constant. The constraint shown would not satisfy Corollary 5.4.1, since it intersects 
dGf inside the cylinder. One way to satisfy the condition in the Corollary is to find the point 
w* of dGf n dGi (if any) for which iuf Pw \ is minimum and enforce wJ*Pw* >1. In other 
words, we wish to solve the optimization problem 


T 

min w i Pw\ 

W2 EM 


[wj I iu 2 ] T € dGf 0 dGi 


s.t. 


(5.18) 


The set dGi can intersect dGf at the half planes TnA w wi = — (TiB w )(±r] — £*) or at w 2 = 0. 
Therefore we solve separate optimization problems for each case and enforce wJ*Pw ] ( > 1 in all 
of them. Noting that the portion of dGf contained in w 2 = 0 is defined by iTji^rci — TiB w (* \ < 

| TiB. w \r] we see that the optimization problems to be solved are 

np 

min w i Pwi s.t. 

«>2 7^0 

TnA w wi = -(TiBuffq sign(G. u ,u>2) - C*) 

Tnwi + V i2 w 2 = 1 


and 


rp 

miri'Wj Pw\ s.t. 
Tjiwi = 1 

TnA w wi - TiB w C\ < \TiB w \rj 


(5.19) 
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Both optimization problems require solving the more general problems 


rp 

minx Px s.t. 
F\x = a 
F 2 x = b 


and 


rp 

minx Px s.t. 
Fx = c 


For the first problem, a unique solution x* is readily obtained by Lagrange multipliers, pro- 
vided that F\ and F 2 are nonzero and non-parallel. The solution is such that x* T Px* = 
[o6]M _1 [a | b] T , where 


M = 


Fi 

F 2 


P~ [F\ 



For the second problem the solution is x* = cF T /(FP~ l F T ) provided F/0. It is straightfor- 
ward to see that M is a symmetric, positive definite two-by-two matrix. Details of the solution 
to the optimization problems are contained in the appendix. Verification of Corollary 5.4.1 is 
achieved through the decision procedure shown in Tables 5.5.4 and 5.5.4. The tables are used 
with the following variables: 


w, 


2,crit 


W. 


2,crit 


Wo 


Wo 


Wi 

Vl 

V2 

V3 

V A 


M22 ~ MnT^Bjr/ - (*) 

r i 2 M 22 

M22 — Mi 2 FjB w (—r) — (*) 

Fj2-Vf22 

1 / | rA(;)-(*) \ 
Ti2 V a ) 

_L( 1 + rjB w (- v -e) \ 

r i2 V a ) 

r nP-^l 

[-TiB w (r,-C) 1] 

[TiB w (r] + C) 1] 
[\TiB w \(v~Q 1 ] 

[-|r.^ 1(77 + 0 1] 


where Mij are the entries of M 1 and Mij are those of M. 


5.5.5 Design Philosophy and Control Constraints 

Unlike [20, 19], The method presented here uses a fixed sliding manifold, which is equivalent 
to fixing a state feedback gain in [19, 20]. That is we do not seek the largest RPI cylinder, as 
this may be achieved with sliding coefficients that do not satisfy other design requirements. In 
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Label 

Condition 

Cl 

to C*) 2 >nfe 

C2 

(”+C) 2 >!T7^ ~ 

C3 

vi M 1 v[ > 1 

C4 

V2M~ Y V2 > 1 

C5 

> M22 

C6 

(r. B .(, + c-))- > 

C7 

TuP- 1 ^ < 1 

C8 

> 1 

C9 

> 1 


Table 5.1: Condition Table 


Case 

Subcase 

Condition 

ft Tji 

Gw w 2,crit > 0 and r *2 / 0 

Cl 

G w w 2,crit < 0 or rj 2 = 0 

C3 

G« J w 2 “ crit < 0 and r i2 / 0 

C2 

G ™ W 2 ,crit > 0 or r i2 = o 

C4 

TiiA^u;! - < TjS w 

C7 

|r a ^u)i - r^C*! > 

C8 and C9 

otherwise 

constraint satisfied 

^ii-A w — ar<i 

G w u >2 > 0 and T* 2 / 0 

C5 

G w W 2 < 0 and T* 2 / 0 

C6 

r i2 = 0 and 7] = C* - a/(TiB w )\ 

C7 

\a - TiB w C\ < \TiB w \r] 

C7 

otherwise 

constraint satisfied 


Table 5.2: Decision Table 
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this case, the method requires that the SMC design be partially completed before addressing 
the constraints. More specifically, the sliding hyperplane needs to be chosen a priori. Then 
the decision procedure is applied to determine an allowable combination of switching gain 
and invariant set. Nominal performance of the sliding motion is thus guaranteed. Several 
refinements can be incorporated to aid the solution of LMI (5.14). For instance, one may rule 
out large ellipsoids, e.g., those whose interior contains the constraint set. This is accomplished 
by enforcing P > kl, where k is a suitable positive scalar to be determined from constraint 
geometry. The volume of the ellipsoid can be maximized under P > kl and the LMI constraint. 
The problem formulation in this case becomes: 


min trace (P) s.t. 
P A w + A^,P + cxP C,' P B\ 

C’BjP -otl 

P > kl 


< 0 


Control constraints are easily incorporated in the design. In fact, the control law of Eq. (5.3) 
can be expressed in x-coordinates as 

u = —GAx — rj sign(s) 


Thus, it is straightforward to show, using the triangle inequality, that a control constraint of 
the form |rt| < u can be accommodated by introducing the additional state constraints T u w < 1 
and — T u w < 1, where 


_ GAJ 

-L u — _ 


U ~ T] 


5.6 Numerical Example 

Consider the following controllable pair: 


' -1 

0 

o 


' 1 ' 

1 

CM 

1 

-1 

II 

cq 

0 

1 

o 

1 

1 

o 


1 


Consider that the constraints in x-coordinates are given by a parallelepiped containing the 
origin in its interior. The rows T x specify individual constraints: 

' 0.8 0.32 -0.16 ' 

-1 -0.4 0.2 

0 0.96 0.32 

0 -1.2 -0.4 

0 0.4 0.8 

0 -0.5 -1 _ 

Suppose that the disturbance is given by ((t) = 0.2 sin(10t). An appropriate set of transfor- 
mation matrices is given by 



' 1 - 

1 

1 ' 


' l 

0 

-1 ' 


' l 

-0.5 

0.5 ' 

T = 

0 

1 

0 

; k = 

0 

1 

-0.5 

; J = 

0 

1 

-0.5 


0 - 

1 

1 


0 

0.5 

1 


0 

-0.5 

1.5 
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The transformed constraints are obtained as the rows of T X J: 


T X J = T = 


0.8 0 0 

-10 0 

0 0.8 0 

0-1 0 

0 0 1 

0 0 -1.25 


Note that the constraints are not symmetric in either x or w coordinates. Using the transfor- 
mation x = TKw = Jw, matrix An has two negative eigenvalues. Choosing to place the poles 
of A w at [—1.5 ± 0.5*] results in G\ = [0 — 0.5], which corresponds to G = [0 0.5 1] in the 
original x-coordinates. The range for a is (0, 3). Choosing an arbitrary fixed a = 2.1, we solve 
the volume optimization problem using k = 0.5 and ?/ + ( = 1. The solution is a matrix P 
with [0.8284, 0.6581] is the diagonal and -0.2278 in the off-diagonal. Choosing r] = 0.8 satisfies 
the decision procedure for constraint qualification and thus the intersection of the cylinder and 
state constraints in re-coordinates is RPI. Of course, the intersection of the transformed cylinder 
and constraints in x-coordinates is also RPI. Using the alternative method with the same P 
matrix, we obtain (3 = 2.3324 and 

Vent, + C = — 7 = = 1.6667 

2 

Then we may choose, for instance, 77 = 1.4 for ellipsoidal invariance alone, but the second 
constraint would be violated. Figure 5.2 sketches the shape of constrained cylinder in w- 
coordinates. Figure 5.3 sketches, in x-coordinates, the trajectories projected onto the X 3 = 0 
plane, and the constraints and transformed cylinder section at X 3 = 0 . 


5.7 Nonlinear Systems 

Certain kinds of nonlinear system under Sliding Mode Control with linear manifolds have 
reaching dynamics similar to that of linear systems. SI nonlinear systems in integrator-cascade 
form are a direct example. Certain minimum-phase, nonlinear MIMO systems under tracking 
control via sliding modes result in decoupled tracking error dynamics which can be also treated 
with the methods described here. In the case of the mixer, such decoupling is possible and the 
output tracking error behaves as 

e p + c p e p = — r/i sign(si) 
e P + c P ep = -r?2 sign(s 2 ) 
e w + c w e w = -7/3 sign(s 3 ) 

where the subindices indicate density, pressure and flow. It is clear that positively invariant 
sets can be found for the tracking errors by considering each equation independently and using 
the methods described above. Valve positions, however, do not appear explicitly in the above 
input/output equations, as they constitute internal states. An analysis of the mixer nonlinear 
dynamics is required to determine if valve positions can be used in a linear state-space realization 
of the above error dynamics. If that were possible, the methods of this work could be directly 
applied. 


56 


RELEASED - Printed documents may be obsolete; validate prior to use. 




Figure 5.2: RPI State constrained-cylinder 


Constraints and Cylinder Projection for x 3 =0 



Figure 5.3: Projected trajectories, cylinder and constraints onto X 3 = 0 
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5.7.1 Nonlinear Systems in Integrator Cascade Form 

Consider the single-input system 


X\ = X2 


X n —1 — % n 

x n = f(x) + g(x)u (5.20) 

where g(x) ^ 0 in a convex region of containing the origin. As in the linear case, let the 
sliding function be defined by s = Gx. Suppose the control input u is such that 

s = —T] sign(s) 


with tj > 0. Assume w.l.o.g. that the n-th component of G is 1. Then the closed-loop dynamics 
can be expressed as 


Xi 


' 0 

1 

0 


Xi 


' 0 ' 

%n— 1 

= 

0 

0 

1 


%n— 1 

— 

0 



_ 0 

—G\ . 

• • -Gn-r . 




_ 1 _ 


= Gx 


7j sign(s) 


(5.21) 


When sliding occurs, x n = — [Gi G 2 ...G n -i].[x 2 ■■■x n ] T and the last state equation in ( 5.21) is 
redundant with the sliding condition s = 0. The reduced dynamics is described by 


Xi 


0 1 

O' 


Xi 

%n— 2 
%n— 1 


0 0 

. —G\ -G 2 . . 

1 

• . -G n _! _ 


%n— 2 
X n —1 


(5.22) 


Define the constant matrix in Eq. (5.22) as A w . This matrix has the standard controllability 
form, and it is clear how to choose the first n — 1 coefficients of G to achieve a stable sliding 
mode. In the following, it will be assumed that A w is Hurwitz. As in the linear case, we seek 
a coordinate transformation which reveals the defective structure of the constant matrix in 
Eq.( 5.21). The existence of such transformation is ascertained in the following 


Lemma 5.7.1. There exists a coordinate transformation x = Jw with J nonsingular in which 
the closed-loop dynamics ( 5.21) are expressed as 



W\ 



0 


w 2 


0 

0 


w 1 
W 2 


+ 


Bi_ 

-1 


r]sign(s) 


{s = w n 


(5.23) 


where B\ = [0 .... 1] T . 
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One possible form of J that verifies Lemma 5.7.1 is given by 


0 1 ... 0 


J- 1 


0 0 
G\ G 2 


1 

1 


(5.24) 


Once the system is written in 7i>-coordinates, it is straightforward to find invariant tunnels with 
the methods of previous sections. 


5.8 Application: Invariant Intervals for Mixer Outputs 

Under the sliding mode control law of Eq.( 4.24), it is clear that the tracking error dynamics is 
given by 


e P + c p e p = -7/i sign(si) 

(5.25) 

e P + c P ep = —r /2 sign(s 2 ) 

(5.26) 

e w + c w e w = —r /3 sign(s 3 ) 

(5.27) 


The developed theory can be directly applied to obtain conditions under which a given interval 
is positively invariant. This is particularly useful in the case of the pressure. If mixer pressure 
is outside the range of supply and discharge pressures, flow reversals and control singularities 
may result. As it will be seen next, enforcing invariant intervals for the pressure will also have 
beneficial effects on tracking performance. Note that absolute temperatures are always positive, 
thus enforcing a positive pressure interval will necessarily result in a positive density state at 
all times and physical feasibility of the control scheme is maintained. 

5.8.1 Invariant Pressure Intervals 

It is straightforward to see that the closed-loop dynamics (5.26) is equivalent to a double 
integrator plant under sliding mode control with manifold 

s 2 = ep + cpe 

Thus the theory presented earlier is immediatly applicable with matrices 

; G = [c p l] 

The objective is to guarantee that P s < P < P\. where it has been assumed, w.l.o.g., that 
Pi =min(Pi, P 2 ). This can be achieved by setting a margin 5 > 0 such that \Pd — P\ = \ep\ < 5 
and by restricting the reference pressure to satisfy Pd < Pi — 5 and Pd > P s + S. Note that this 
is in fact a performance requirement, since bounds are being placed on the tracking error \ep\. 
Although the theory so far developed will yield invariant sets of simple description in some 


A = 


0 1 
0 0 


B = 


59 


RELEASED - Printed documents may be obsolete; validate prior to use. 



set of coordinates, it turns out that, in this case, positively invariant intervals in the original 
coordinates are also directly obtained. The transformations T and K and J are chosen as 


T = 


1 0 
1 1 



The reaching dynamics is expressed in 


1 -1 

-(cp + 1) 1 

re-coordinates as 


; J = 


l -l 

— cp 0 


w\ = — cpW\ H rj sign(s) 

cp 

W r 2 = — rj sign(s) 

c P 

s = — 10rc2 


( 5 . 28 ) 


Since \e\ = cp\w\\, it can be readily established that the interval defined by 


invariant if 


V — 


c 


3 

p 


5 


< 6 is positively 
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Chapter 6 

A Graphical User Interface for 
Mixer Simulations 


The fairly complete modeling and control analysis effort for the mixer must be complemented 
with a simulation environment. The purpose of this activity is to develop and program a Graph- 
ical User Interface (GUI) for mixer simulations using the validated model and its controls. The 
objective is to generate a GUI which allows simulations and control tuning without advanced 
knowledge in control theory or mathematical models, so that operations personnel can benefit 
from the results of this research. While the programming details are left out for not constituting 
a research activity, some of the design criteria and a basic functional description of the GUI 
follow. Although some packages like EASY5 and GFSSP [14] are available that incorporate the 
basic laws of fluid flow and thermodynamics, the objective here is to construct a simulation 
environment that trades off generality by increased flexibility, modularity and ease of use in 
test stands and associated technologies. Key innovations are cited below 

1. Use of advanced concepts in modeling and control: Test stands and related subsystems are 
fairly predictable and governed by basic laws of thermodynamics and fluid flow. Therefore 
reliable dynamic models are possible to construct. Models generally consist of systems 
of nonlinear ordinary differential equations. Model properties that are relevant to ac- 
tual system operation can be inferred through mathematical analysis of the model. The 
determination and tuning of control strategies can be likewise carried out with the aid 
of modern control theories. The currently available hardware does not permit imple- 
mentation of advanced concepts in real-time operations. That is, most real-time control 
processing is implemented in programmable logic controllers (PLC). Advances in process- 
ing speed and hardware/software architectures by 2006 will be suitable for the deployment 
of advanced concepts in modeling and control. 

2. Modularity and expandability: Isolated subsystems are modeled in a modular fashion. 
Each model has a number of inputs and outputs and corresponding connection rules. 
Modularity allows the simulation of individual subsystems by feeding them with pre- 
defined data, or of larger systems by interconnecting the inputs and outputs of various 
subsystems according to connection rules. Modules may be added indefinitely to create 
expanded systems. An important advantage of a modular set of models is that a working 
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simulation is available as soon as the first module is complete, thus making possible a 
concurrent use and development of the simulation tool. 

3. Flexibility and user-friendliness: Programs such as EASY5 and GFSSP [14] can, in prin- 
ciple, be customized to meet simulation needs for test stand facilities. They, however, 
lack the flexibility that the proposed environments will provide. Both packages require 
a trained and experienced user. Moreover, including feedback loops may require the 
writing of code. The proposed simulation environments will present graphical user inter- 
faces (GUI) that mask the mathematics governing the model in a black-box approach. A 
user needs to be familiar with the actual test stand only, not with a particular scripting 
language. 

4. Accuracy: The incorporation of the latest modeling and control technologies in the mod- 
els and simulation environments will produce more accurate simulations. Future test 
requirements will demand higher accuracies than those possible with current prediction 
and control schemes. The modeling work proposed here will satisfy future demands for 
higher prediction accuracy. Accurate simulations and model validation has a direct im- 
pact in design and operation of test stands, and both activities will need to be done more 
effectively than today. 

5. Integration: The proposed simulation tools will link with health management technolo- 
gies (HMT) directly. The simulations will provide inputs to HMT and get important 
information about system changes from HMT. 

6.1 Functional Description 

The front interface is shown in Fig. 6.1. A brief description of GUI functions follows. 

• Input Fluid Properties: The user enters either a constant pressure and temperature for 
the GH2 and LH2 or clicks on the From File button to import two ASCII-formatted data 
files. 

• Mixer Parameters: The simulation must start with some pressure and temperature of the 
fluid inside the mixer. These constants, along with the volume of the mixer in cubic feet 
are entered. 

• Exit Pressure: The model considers the pressure past the exit valve a boundary condition, 
that is, it must be supplied as a parameter. The user has the option of specifying a 
constant pressure by directly typing it in the box, or using a time-varying pressure from 
a file. 

• Control Valves: Each control valve is defined in the same way. The user clicks on the 
Parameters button and enters the value of the flow coefficient when the valve is fully 
open, the time constant, and the gain between percent opening and commanded opening 
(in percent per volt). In this version, the relationships between flow coefficient, percent 
opening and commanded opening are taken as linear. The user then clicks on the Source 
button and specifies how the particular valve will be operated. The choices are to leave the 


62 


RELEASED - Printed documents may be obsolete; validate prior to use. 



File About... PID Controls Calculators Operating Sequences.. 


X 

□ 









m 











Figure 6.1: The front interface 
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percent opening fixed by typing a value between 0 and 100 in the Constant box, to specify 
a ramp by entering the appropriate values, or to import a valve opening history from a 
file. The file must be in ASCII format and contain two columns of data, the first one 
containing time and the second one containing the corresponding percentage opening. If 
the Enable Dynamics box is checked, a first order transfer function with the specified gain 
and time constant is assumed between commanded and actual valve openings. Otherwise, 
the GUI takes a zero time constant and sets the gain to one. 

• Profile Tracking: This area is used to enter desired output profiles when operating un- 
der feedback control. The user chooses the outputs to be constant setpoints, ramps, or 
arbitrary profiles imported from a data file. 

• Control Strategies: The user selects a control option from a pull-down list. The following 
options are available: 

1. Pressure PI Control (GH2 Valve Only): This simulates the automatic stage of the 
control mode often used at the test stands. A feedback loop is formed by measuring 
mixer pressure and using the tracking error in a proportional-integral law. The GH2 
valve is actuated, while the LH2 and exit valves remain at fixed positions. The 
proportional and integral gains are selected in the PI Controls pull-down menu at 
the top-level interface. 

2. Multi Input PI Control: This control law simultaneously actuates the three valves. 
It is based on the work done by Barbieri [3]. The control law was derived based on 
local linearization of the model equations. It is useful for regulation close to a steady 
setpoint or for commanding the output variables in small ranges. 

3. Feedback Linearization: This control law simultaneously actuates the three valves. 
It is based on the work done by Richter [26]. The control law may be used for 
regulation or tracking over wide ranges of the output variables. The gains involved 
in the control law are hidden from the user in this version. 

4. Manual Valve Control: If this option is selected, simulation proceeds in open-loop, 
using the valve profiles specified in the Control Valves area 

5. Operating Sequence: This option allows to schedule switching between control strate- 
gies during a simulation. Simulation may be initiated, for instance, in manual mode. 
When switching criteria are met, the control mode is switched to one of the auto- 
matic options. The automatic mode may be then released and operation may be 
switched to a shutdown procedure. Switching criteria is entered as time values or 
values for combinations of process variables. 

6.1.1 Calculator Functions 

In addition to the dynamic simulation modes, the GUI operates as a steady-state calculator, 
namely two functions are incorporated: 

• Setpoint Targeting: The GUI calculates the constant valve openings required to obtain 
the steady output values specified as Constant in the Profile Tracking area. The results 
are displayed in the text boxes. 
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• Operating Point: This operation is the inverse of Setpoint Targeting. It calculates the 
steady mixer pressure, exit flow and exit temperature that result when the control valves 
are left at the constant openings specified in the Control Valves area. 

6.1.2 Operating Sequences 

Actual mixer operation is frequently done according to a schedule where switching between 
manual and closed-loop controls may occur. In one of the typical operating modes, the mixer 
is first manually brought close to a desired operating point and a feedback controller is then 
engaged. The Operating Sequence feature allows to simulate such modes. The user specifies 
control strategies, control objectives and manual valve actions using the interface. Upon select- 
ing the Operating Sequence mode, the user is prompted for a condition to be satisfied for the 
current control settings to be effective. The user enters the condition as a Matlab string, for 
example ’ (t>2) & (abs (P-7000) ) <=10 ’ specifying that the current control settings become 
active if the pressure is between 6990 and 7010 psia and the time is greater than 2 seconds. 
The GUI then records the settings and adds them to a schedule. An arbitrary number of events 
can be scheduled. If Operating Sequence is selected as a control strategy, then the simulation 
proceeds according to schedule. 

6.1.3 Simulation Control and Results 

In this version, the user has control over the simulation time span, and can start and stop a 
simulation by clicking on the appropriate button. The simulation is performed by Simulink. 
Advanced users may access the file containing the simulation diagram and modify parameters 
such as numerical method, time steps and tolerances. If the Plot check box is enabled, the 
simulation results are displayed in windows, from where they may be printed. 
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Chapter 7 

Pressurization and Collapse in 
Cryogenic Tanks 


The high oxidizer and fuel pressures typically required by engine and component testing are 
achieved by means of a run-tank equipped with a pressurization mechanism. A typical spherical 
run tank is depicted in Figure 7. The tank is partially filled with the liquid propellant to be 
pressurized. The space above the liquid surface is occupied by a gaseous pressurant, typically 
Nitrogen for LOX run tanks and Hydrogen for LH2 run tanks. The gas contained in the tank is 
known as ullage. There exist several methods of pressurization [28], broadly classified into self- 
pressurization and external pressurization. Self-pressurizing systems use the vapor produced 
by the boil-off of the liquid propellant due to heat transfer. In this project we concentrate 
on external pressurization systems, which rely on an independent source of pressurant gas. 
In this configuration, the liquid is pressurized by controlled injection of pressurant from high 
pressure bottles. A control valve is installed between the upper portion of the tank and the 
storage bottles. This valve serves as primary control element for the tank pressurization. The 
pressure is usually controlled with a PI loop as shown in Figure 7.2. Pressure collapse of cryogen 
propellant during pressurization has been identified as an effect that often leads to down time 
and test disruptions. Physically, the collapse occurs when the incoming stream of warm gas 
is violently cooled by contact with the cryogenic liquid. The sudden drop in temperature 
creates a pressure drop which counteracts pressure build-up from the added mass. Depending 
on key parameters such as incoming flow rate and temperature and degree of mixing, the two 
competing effects result in an initial pressure build-up followed by collapse. In some cases, 
pressure recovery is observed. The phenomenon is worsened by control valve saturation and 
the subsequent integrator windup in the PI loop. The purpose of this section is to construct the 
simplest possible dynamic model, consisting of ordinary differential equations, that captures the 
collapse phenomenon qualitatively and with an acceptable degree of numerical accuracy against 
activation data. Once the model is validated and refined, it is desired to study the effects of 
the control loop on the collapse phenomenon, initially by simulation with various PI gains 
(including actual gains) followed by an in-depth analysis to determine if a properly designed 
controller can mitigate the collapse effect. 
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Figure 7.1: Oxygen run tank schematic 


Pressure Setpoint 



^ Pressurized LOX 
to test article 


Figure 7.2: Run-Tank Pressurization System 
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7.1 Simplest Model for Tank Pressurization 


A fairly detailed model of the run-tank that feeds the mixer was developed by Follet and 
Taylor [7] and programmed in EASY5/ROCETS for simulation studies. While the model is 
adequate to predict the behavior of variables related to tank operation, it introduces complexity 
that makes it difficult to use it for controller design purposes. Specifically, the calculation of 
heat transfer through tank walls and between the vapor-liquid interface is what introduces 
complexity in the form of cumbersome mathematical expressions and iterative procedures. As 
stated by the authors in page 6 of [7] , the heat transfer correlations use Rayleigh numbers that 
exceed the recommended limiting values, possibly introducing errors. In the following model, 
the heat transfer through tank walls is set to zero, while the exchange of heat between the 
phases is not explicitly considered in the equations, for the energy balance equation uses the 
whole tank as a control volume. Mass conservation is still applied separately to the liquid and 
vapor phases, and the assumption that there is no net exchange of mass through the vapor- 
liquid interface is also used. The model in [7] uses the same assumption with good results. The 
proposed control model will be compared to the existing model in its prediction capabilities. 
The basic governing equations can be obtained from conservation of mass and the First Law of 
Thermodynamics for a control volume with transient terms. Refer to Figure 7.2 for a depiction 
of the pressurization system. The following assumptions are made in the derivation of the first 
model: 

• The tank is adiabatic. 

• A single coefficient is used to describe the heat transfer between gas and liquid under all 
conditions of stirring and mixing. The coefficient, however may be allowed to vary during 
pressurization. 

• No gas is condensed. 

• The inlet valve behaves isoenthalpically. 

• The ullage and liquid pressures are equal at all times. 

• The bulk temperature of the liquid does not change during pressurization. 

Denote the source gas mass flow rate and enthalpy by w and h\. respectively. Let T, p. P and 
u denote the ullage gas temperature, density, pressure and internal energy. Let T c be the bulk 
temperature of liquid and V (t) denote the instantaneous volume the liquid in the tank, so that 
the work done by the ullage is given by V(t)P. The energy and mass balances for the ullage 
give the following system of differential equations: 

P = {«>(*) - T(t)p} (7.1) 

u = [~V(t)P + w(t) (hi -u)- r](T - T c ) } (7.2) 

where r] is a heat transfer coefficient which encompasses all mixing effects. Note that all 
quantities -except T c - are time- varying, however the explicit time dependence notation has been 
used for V (t) and w(t) only, to emphasize that they are exogenous inputs to the dynamic system. 
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According to [28, 15], p can be calculated from a natural convection correlation involving the 
Prandtl and Grashoff numbers of the ullage gas: 

r] = h c A, where: (7.3) 

h c = K h C^{ GrPr) n 

Kh, C and n are constants, Kf is the thermal conductivity of the ullage, d is a characteristic 
length, A is the area of the liquid-ullage interface, and Gr and Pr are the Grashoff and Prandtl 
numbers, respectively. The above equations already show a basic qualitative description of the 
competing effects: as more mass enters the tank with rate w, the internal energy increases due 
to the second summand in Eq. (7.2). However, the third summand is negative whenever the 
ullage gas is warmer than the liquid, creating a counteracting effect. The interplay of these 
contributions to the energy derivative depends on the operating parameters, mainly rj, T c and 
w. Although the description given in [28] is mostly academic, it is possible to leave rj as an 
unknown function of to be identified from experimental data. The identification of p is listed 
as future work in Section 8.1. Note also that the pressure P and temperature T of the ullage 
are functions of the state variables p and u, as are the variables participating in the calculation 
of 7j and w(t), making Eqs. (7.1) and (7.2) a self-contained system of differential equations. 
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Chapter 8 

Conclusions and Future Directions 
for Research 


Tenure lasted 31 months, during which a significant amount of results relevant to rocket test 
stand modeling, simulation and control have been produced. Moreover, the research activity 
has inspired novel theoretical work in the area of constrained control and has revealed new 
directions for research at the intersection of control theory, thermodynamics and fluid flow. 

8.1 Work in Progress and Future Work 

Modeling work on test stand subsystems can benefit from the development of an efficient ther- 
modynamic property library for Matlab, relying on correlations instead of table look-ups and 
interpolation. Such work has been completed for Oxygen and will be continued for Hydrogen 
and Nitrogen. 

The theoretical developments in constrained sliding mode control need to be expanded to 
include the multivariable case, to consider tracking problems and to incorporate a chattering 
reduction mechanism which is compatible to the cylinder methodology. The boundary layer 
approach of Slotine [30] seems appropriate, as the boundary layer itself enjoys invariance prop- 
erties. The tracking errors in the mixer problem have linear dynamics falling in the framework 
of the cylinder approach, however the possibility of applying the theory to constraints in the C v 
coefficients is still unclear. That is, the theory can be used to limit the ranges of thermodynamic 
variables, but the fundamental issue of valve position constraints is yet to be solved. 

Work on the run tank model requires experimental data tailored to the needs of the project. 
It remains unclear at this point if the data routinely collected during tests is adequate for 
system identification and parameter estimation procedures. Upon a successful validation of the 
run tank model, its interconnection with the mixer can be modeled. 

Work on the GUI can continue by developing an interface for the validated run tank model. 
Ultimately, a collection of GUIs developed for independent subsystems can be brought together 
in a GUI suite which facilitates large-scale simulations. While the GUI presented here is proto- 
typical, future versions need to address the issue of execution speed, which is already low, even 
for the mixer alone. Compilation of modules and routines into dynamic run-time libraries may 
improve execution speed. Substitution of the look-up routines for thermodynamic properties by 
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efficient ones based on correlations is also expected to significantly improve execution speed. 
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Appendix A 


A.l Construction of T 


One of many ways to select T is presented here. Partition T congruently to A and B and write 
the requirements as 

T\\b\ + T1262 = 0 
T ‘ 2 \ b \ + T22&2 = 1 

where the T t] denote the block components of T^ 1 and b\ and 62 those of B. Denote the 
( 1 , 1 ) component of T~ 1 AT by An, as in Section 5.3 and let a %] denote the components of A. 
Suppose, w.l.o.g., that the scalar 62 is nonzero, which can be always achieved by re-labeling the 
states if necessary. Then select any invertible Tn and choose 

rr Tubi 

Tn = -~ 


Selection of T21 can be done so that the ( 1 , 1 ) block of T~ l AT has desired eigenvalues, specif- 
ically, to make it nonsingular. In fact, using the block matrix inversion results from [ 11 ] and 
performing lengthy algebra it is possible to show that An = TnAnTn \ where An = A — BT21, 
with 


A — an~ 


M21 

b 2 


and 


B = 



bi + 



b 2 


Thus T-21 can be used to assign nonzero eigenvalues to An with the restriction T‘2\b\ / 1. Since 
An is similar to An, it is nonsingular. The selection of T is finished by choosing 


T22 


1 — T 2 1 b 1 

h 


Note that T22 7^ 0 , therefore T is indeed invertible [ 11 ]. 
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A. 2 Proof of Lemma 5.3.1 


Since 

(I - BG)A = T 

J can be chosen as TK where K is an invertible matrix such that 



Using the block matrix inversion results of [11] it can be shown that the following is a solution: 

K = \ 1 zAuvlu 

L"~— Gi i 

Note that invertibility of An is guaranteed by choice of T as explained in Section A.l. More- 
over, with such choice of J we have that the sliding manifold coefficients are expressed in 
^-coordinates as 

s = GJw = (GT)Kw = [Gi | 1] Kw = [0 | - G\A^Ai 2 + l] w = G w w 2 
Also, B 2 can be readily found to be 

B 2 = (G^Au - l)" 1 

which results in B 2 G W = — 1. 

A. 3 Proof of Theorem 5 

(Sufficiency). In view of Th. 2, we wish to show that f(w, () £ Cs(u>) Vrc £ <95, VC £ Z. It 
is possible to show, under the convexity and closure assumptions, that one can decompose the 
boundary of S as 

dS = (M n 8G) U (dM n G) 

First, we show that f(w, C) £ Cs(w) V re £ ( M n dG ). By hypothesis, Tif(w, Q < 0 Vic £5(1 
dGi, \/( £ Z, i = 1, 2 ..m. Noting that |J(Sfl<9Gy = MC\8G we have that f(w, C) £ Cg(w)Vw £ 
(M n dG). Now, since 8G fl M C M we have that, if w £ int(Af) then Cm(w) = thus 
f(w,C) £ Cm(w), whereas if re £ 8M a then f(w,Q £ Cm(w) due to M being RPI. Therefore 
f(w, C) £ Cm(w) FI Cq(w) M8G n M. By assumption, 0 £ int(M — G). Then, by Property 5.2.1, 
Cm(w) FI Cg(w) = Cmdg(w) = Cs(w). Thus f(w,() £ Cs(w) \/w £ ( M fl 8G). Now we show 
that f(w, C) £ Cs Vw £ (Gfl 8M ) . M is RPI by assumption, therefore f(w, C) £ Cm Vw £ 8M , 
in particular Vw £ 8M fl G. Since 8M fl G C G we have that, if w £ int(G) then = M n , 
thus f(w,Q £ Cg(w), whereas if w £ 8Gi for some i then w £ 8M n 8Gi C 8Gi fl M, since 
M is closed. By hypothesis, it follows, as before, that f(wffi) £ Cg(w). Thus f(wffi) £ 
Cm(w) Fl Cg(w) = C$(w) \/w £ 8M fl G. We conclude that f(w, () £ Cs \/w £ 8S , VC £ Z , and 
that S is therefore RPI. 

(Necessity) . Suppose S is RPI, and by contradiction, suppose 3vj £ S fl 8Gi , 3C £ Z for some 
i such that Tif(w, () > 0. Then f(w , C) ^ Cg^w) fl Cm(w). So 3w £ S fl 8G , 3( £ Z such that 
f(w,() ^ Cq(w) fl Cm(w) = Cs(w). Since S fl 8G C 8S we have that 3w £ 8S , 3( G Z such 
that f(w, C) ^ Cs(vj), contradicting that S is RPI. 
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A. 4 Proof of Lemma 5.4.1 


Partition P,; as Fj = [P,;i | P^] • Let dGi = {[ref | W 2 ] T £ M n : Tnw\ + r ? ;2rc2 = l}. Suppose 
T i2 / 0. Then dGi can be parameterized as j[ref 1 1 ^ 2 ] T : re 2 = — , re 1 £ M n_1 j. Con- 
sider the set V = jre 2 £ M : W 2 = — , rei £ <S 0 j. The functional h : M n_1 — > M given by 
h(wi) = ; yM 1 rci is bounded and therefore continuous in M n_1 . Then h(S a ) is compact, since 
S 0 C M n_1 is compact [16]. Thus V, the translation of h(S 0 ) is also closed and bounded, i.e., 

V is a real closed interval posessing a minimum and a maximum. Note also that the maximum 
and minimum values of V are achieved at the boundary of S Q . Denote by ref and the 
half-spaces of M n where W 2 > 0 and W 2 < 0, respectively. Denote by re 2 the plane W 2 = 0. 

Consider, first, the case when dGi PI TLs 0 [w 2 , ^ 2 ] Pi w 2 = 0. Let b =min(F) > 0. Note that 
w 2 > b > w 2 . Then 11.sJ(UL 2 ,b\ is RPI as a member of the family TLs 0 - This implies, by Th. 2 
that f(w, C) £ CH So \w 2 ,b}( w ) Vw € &Hs 0 \w 2 ,b], V( £ Z. Now, Hs 0 {w 2 ,b] C 2 ], 

thus C Hso \^ 2 ^(w) C C GnWso [^ 2i5j2 ](u;) in particular for a point w b = [ref \ b] £ dHs 0 [w 2 ,b], 
where w x £ dS a . So f{w b ,() € C GnHSo[m ^ 2] (w b ) = C Gi (w b ) nC HSofeif2 |(ro t ). Therefore 
L ? :/(rcfe,C) < 0, with w b £ dGi Pi dHs 0 \w 2 , W 2 ] C PI dis 0 [w 2 ,u) 2 \- To prove the strict in- 
equality, suppose, by contradiction, that Fif(w b ,() = 0 for some ( £ Z. This would require 
r i2 B 2 (±v — C) = 0) which is impossible, since £ < ( < ij. Thus w b is a point satisfying the 
Lemma. The case when dGi Pi B.s 0 [w 2 ,W 2 ] Pi ref = 0 is treated similarly, taking b =max(V), 
which is negative, and considering the cylinder 'Hs„ [b, w 2 ] ■ Finally, consider the remaining pos- 
sibilities of dGi Pi H So [m 2 , ^ 2 ] Pi w 2 / 0 or V i2 = 0. In those cases, the intersection of dGi 
with the cylinder and the sliding hyperplane W 2 = 0 is nonempty. System dynamics on the 
sliding hyperplane is given by w\ = A w w\. Consider the set resulting from the intersection 
of S 0 and the closure of the complement of the state constraint restricted to w 2 = 0, that 
is, let K = S 0 Pi G c oi , where G c oi = {w\ £ M n_1 : T^rci > l}. It is easy to see that K is 
compact and convex and that it does not contain the origin. Suppose, by contradiction, that 
Turhi > 0 V w\ £ dG 0 i Pi S a . Then, following arguments similar to those used in the proof of 
Th. 5, it is possible to deduce that K is positively invariant for w\ = A w w 1 , which contradicts 
the asymptotic stability of the origin. Thus Elrei £ dGi a Pi S G such that T 0 iW\ < 0, which in 
turns implies that re = [ref | 0] £ dGi Pi Bs 0 satisfies r^rc < 0. 

A. 5 Proof of Proposition 5.4.1 

For each re £ M n , consider the set a(w) = {Tif(w, () : ( £ Z}. We show that max a (re) = 

Tif(w,(*), where C* =-sign.(FiB w )( and thus, Tif(w,(*) > Tj/( re,£) \/( £ Z. First note that 
r,:/(re, £*) £ a (re) since C* £ Z. For any ( £ Z we have that 

L ? :/(re,C) = Fi 1 A w wi+F il B 1 (r]sign(G w W 2 )-()+Fi 2 B 2 (r]sign(G w W 2 )-() = F il A w w 1 +F i B w (r]sign(G w w 2 )-() 
and thus 

Tif(w,C) -F ? :/(re,C) = \FiB w \( + F l B w ( 

But since |£| < (we see that F i.B w ( < |Fj.B w ||<f = ir^i^CI < iTt-B^IC which shows that 
|ri.B w |£ + FiB w ( > 0 and therefore max a(re) = Fif(w,(*). To show that the collection of 
subsets is a partition we must show that the subsets are not empty, disjoint and that their 
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union equals R n . Assume, w.l.o.g., that G w > 0. Take w = [ref | W2 f r in the half-space W 2 > 0. 
Since TnA w / Owe have that Tif(w, £*) = 0 is a half-hyperplane separating the half-space into 
open regions for which either T if(w, £*) > 0 or Tif(w, (*) < 0. Any w contained in the positive 
region of the half-space belongs to Gf, thus int(GU) / 0. Similarly, any w contained in the 
negative side satisfies Yif{w,Q*) =maxa(w) < 0. Then Tif(w,Q < 0V£ G Z and thus w G G ~ 
and int(GT) ^ 0. Also, any point w on the half- hyperplane is an adherent point of int(Crf) 
that is not contained in int(GU), thus w G dGf / 0. If we take the half-space W 2 < 0 we can 
likewise show that any point w on the corresponding half- hyperplane satisfies w G dGf. We 
have shown that dGf D { [wf | W 2 ] r G R n : = 0 , W 2 / 0}. Now suppose W 2 = 0. 

We show that dGf D {[tcf | 0] T G R n : \TnA w w\ — TiB w (*\ < \TiB w \r]f Let w\ G R" such 
that \TnA w w\ — TiB w (*\ < \TiB w \r] and e > 0 and consider the ball B([wf | 0] T ),e). Choose 
V 2 such that |T,; | < rjfi u ,sign(G^r; 2 ) and |u 2 | < £, which can always be found. Consider the 
point v = [wf | V 2 ] T ■ Using the euclidean norm in R n we see that ||r> — [wj 1 | 0] T || = |u 2 | < £, 
thus v G B{[wf | 0] T ),e). Moreover, 


TnA w wi -TiB w C* > -\TiB w \r] > - TiB w r)sign(G w V2 ) 

which leads to Tj/(u, C*) > 0, so v G Gf and w = [wf \ 0] T is an adherent point of Gf. However, 
w ^int {Gf), for all the balls B(\w[ \ 0] r ),e) also contain the point v' = [wf \ — V 2 ] T , for which 
C*) < 0. This implies w G dGf, showing the desired inclusion. We have thus shown 

that 

dGf D { [wf I w 2 ] T G R n : Tif(w, C) = 0,w 2 ^ o}u{ [ref | 0] T G R n : paA^ - T Z B W C \ < \r t B w 

Now we prove the reverse inclusion. Suppose that w G dGf and that Tif(w,C*) < 0. Then 
Tj f(w, C) < 0 \/( G Z and w G Gj . Observing from the definition that Gf 0 G~ = 0, the only 
possibility is that w G int(Gf), which contradicts the hypothesis. Therefore Tj/(ie,£*) > 0. If 
Tif(w,C*) = 0 then we have proved the inclusion. The last possibility is that > 0. 

If we take w in either half-space W 2 > 0 or W 2 < 0 we see that the corresponding half- 
hyperplane separates the half-space into two open regions. Then w cannot be a boundary point 
for Gf. Thus it must be that W 2 = 0. By contradiction, suppose \TnA w w\ — TjH !i ,C’ , '| > 

\TiB w \rj. Then take the open ball centered at w that contains all points v = [vf \ V 2 ] T 
for which \TnA w v\ — TiB w C*\ > \TiB w \r). Take v in the ball and suppose V 2 / 0. Then 
\rj(v,C)\ = \TnA w vi + TiB w (±7] - C*)\ > \\TnA w vi - TiB w C\ ~ \TiB w \r]\ > 0. When 
V 2 = 0 we have that |Tj/(u,C*)| = \Yi\A w v\ — | > |TjH w | > 0. Since > 0 

and iu is in the ball, we see that Tif(v,(*) > 0 and thus the ball is entirely contained in 
Gf. This contradicts w being a boundary point of Gf. Thus we have shown that dGf = 

{[wf | W2] t G R n : Tif(w,C) = o , W2 / 0 }u{[wf | 0] T G R n : |T a A w u;i - TiB w (*\ < \TiB w \rj). 

We are now in a position to show that the sets are disjoint. From their definition, it is obvious 
that mt(Gf )nint(G“) = 0. Also, dGff]mt(Gf) = 0. We now show that dGf nint(GT) = 0. 

Suppose w G dGf. Then w /int(GT), since every ball centered at w must contain points of 
int(Gf ) so there is no ball entirely contained in int(G~). It is now straightforward to show that 
int(G+)Uint(G“) U dGf D R n . 
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Appendix B 


B.l Solution to Optimization Problems 

B.1.1 Calculation of f3 (Section 5.5.1) 

The problem is 

minx T Qx s.t. 

X Px = 1 

Define the augmented Lagrangian function as 

C = x t Qx + A (x T Px — 1) 

The first-order conditions are 

dC 

— — = 2 Qx + 2 A Px = 0 
ox 

dC rp 

— — = x 1 Px — 1 = 0 
o A 

The first condition can be written as 


{P~ l Q + XI) x* = 0 


after multiplication by P _1 . Thus —A and x* must be an eigenvalue-eigenvector pair for P^ 1 Q. 
A feasible solution set of vectors x* is obtained by scaling the eigenvectors of P^ 1 Q so that 
they “fit” in the ellipsoid defined by P. If x* /'y is used, with 7 6 M, then the scaling factor 7 
must satisfy 

x Px = 7 

All eigenvectors x* of -P -1 Q are found, and /3 is chosen as 


/ 3 = min 


x* t Qx* 

X*T JP 
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B.1.2 Enforcement of Corollary 5.4.1 


It is straightforward to note that if T *1^4^ = 0 or YiB w = 0 then the constraint is automatically 
satisfied. Assume fji Au, ^ 0 and YiB w ^ 0 and define the following variables: 


M 


Ei 


Fi 

F 2 


P-'[Fj 


Fn A 


W 



F 2 = Tn 


Problem 1: w 2 =/=■ 0 

First, suppose YnA w ft r*i and r* 2 / 0. The problem to be solved is 

min w j Puj i s.t. 

F\W\ = a 
F 2 w\ = b 


with a = — (TiB w )(r] sign(G w w 2 ) — £*) and b = 1 — Yi 2 w 2 . The solution is readily obtained by 
Lagrange multipliers and results in a unique minimizing argument w\ such that 

wf r Pwl = [a b] M _1 [a b] 1 = fi(w 2 ) 


Suppose G w w 2 > 0. Function f\ is clearly quadratic in w 2 and posesses a unique global 
minimum. Differentiating and equating to zero gives 


dfi(w 2 ) 

dw 2 


2 [0 - r i2 ] M _1 


-r^fa-c*) 

i - r i 2 w 2 


= o 


This results in the linear equation 

Y 2 2 M 22 w 2 = T i2 (M 22 - M 12 TiB w (ri - (*)) 


where M, t j are the entries of Af 1 . Note that M 22 ^ 0 since M is symmetric and positive- 
definite. Let W 2 Crit be the solution to the above equation: 


+ _ M 22 — Mi 2 TjB w (r] — (*) 

W2 ' crit ~ r i 2 Af 22 

If G w W 2 Crit > 0 and W 2 crit € (w 2 , tc 2 ), the solution is valid and we enforce vM~ 1 v T > 1, where 
v = YiB w (iq — £*)[—! J^\ ■ After some algebra, this reduces to 


(v-C) 2 > 


M u 

{YiB w y 


If G yj w 2 c rit T < 0 and w 2crit £ (w 2 , w 2 ), the minimum is obtained at G w w 2 = 0 + and the 
condition to enforce becomes vM~ 1 v T > 1, with v = [—YiB w (rj — £*) !]• An entirely analogous 
procedure is followed for the case G w w 2 < 0. Let 

_ _ M 22 - M 12 YiB w (-rj - C) 

W2 ’ crit Y i2 M 22 
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If G w w 2 crit < 0 and w 2 crit G (w 2 , w 2 ), we enforce 


(n + C) 2 > 11 

w ; j (r^) 2 

If G w w 2)C rit— > 0 and w 2crit G {w 2 ,w 2 ), the condition to enforce becomes vM~ l v T > 1, with 
v = [TiB w (r] + £*) 1]. Now suppose r^yl^ F,; i but = 0. Once a sign for G w w 2 has 
been chosen, the function f± becomes constant. Therefore the condition to enforce becomes 
vM~ 1 v t > 1 for both v = [TiB w (r] + C*) 1] and v = [—TiB w {rj — £*) 1]. Now suppose 
TjiAu; = qF,;i and Y l2 / 0. Note that a < 0 is the only possibility for A w Hurwitz. In this 
case, we reduce the restrictions as follows. Substitution gives the equation 

a(l - T i2 w 2 ) = -TiB w (r] sign {G w w 2 ) - £*) 

First we seek a solution for G w w 2 > 0. Define 

+ W 1 + rAh-n\ 


If G w w 2 > 0 and w 2 G (w_ 2 ,w 2 ), then the solution is valid and the reduced optimization 
problem becomes 


F 2 w\ = - 


min w j Piu i s.t. 

^(,-n 


The unique solution w * is readily obtained by Lagrange multipliers and is given by 

r 


The condition to be enforced reduces to 


TiB w (ri-C) 


> M 22 


If G w w 2 < 0 or w 2 (w 2 ,w 2 ) the solution is not valid and the constraints have empty 
intersection. An entirely analogous procedure is followed for the case G w w 2 < 0. Let 


If G w w 2 < 0 and u> 2 G (w 2 ,w 2 ) then the solution is valid and the condition to be enforced 
reduces to 


YiB w (r] + £* 


> M 22 


If G w w 2 > 0 or if w 2 ^ ( w 2 ,w 2 ) then the solution is not valid and the constraints have empty 
intersection. Now suppose Tji A w = aF, i and Ti 2 = 0. Then substitution gives the equation 
a = — TiB w (r] sign(G u ,w 2 ) — £*). If r] = |£* — a/(YiB w )\ then the constraints reduce to Fauq = 1 
and the problem becomes 

ruin w l Pw\ s.t. 

F 2 w\ = 1 

The solution is straightforward and reduces to the enforcement of rjiP _1 r? 1 < 1. 
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Problem 2: w 2 = 0 

The problem is 


min wj Pwi s.t. 

Tnwi = 1 

\TnA w wi - T;.B W C*| < \TiBw\ 

The inequality can be written as 

-| TiB w \{j] + C) < TnA w wi < | TiB w \(ri - () (B.l) 

Note that the interval specified for TnA w w\ is nonempty, since 77 > 0, ( > 0 and 77 > (. The 
problem is of the form 

rr ’ 

minw 1 Pwi s.t. 

F 2 w\ = 1 
b < F\W\ < a 

with a > b, F 2 / 0. The associated problem 

min wj Pw± s.t. 

F 2 wi = 1 (B.2) 

posesses a global and unique solution tcJ. If the solution satisfies inequality (B.l) then w\ is a 
valid solution and must be used in the enforcement of w^ T Pw* > 1. If not, the following two 
problems must be solved 

rp 

min w j Pw \ s.t. 

F 2 w 1 = 1 

F\ wi = a (B.3) 

and 

min wj 1 Pw \ s.t. 

F 2 wi = 1 

F\w± = b (B-4) 

and wl T Pwl > 1 must be satisfied by both solutions. Regardless of whether TnA w || T,i, the 
solution to problem (B.2) is given by 

u .*_ P ~ lT « 

1 TnP-iTl 

If w\ satisfies inequality (B.l), then the condition to be enforced is T < 1. If not, we 
solve problems (B.3) and (B.4) using a = \TiB w \(ri — C*) and b = —\FiB w \(r] + C*)- Suppose 
first that TnA w T,; 1 . Then the solutions to problems (B.3) and (B.4) result in the condition 
to be enforced: vM~ 1 v T > 1, for both v = [\TiB w \(r] — () 1] and v = [— |T,;i ?^|(?7 + () 1], Now 
suppose that TnA w = aT, 1 . Substitution of this into the restrictions results in the equations 
a = | TiB w \(i) — C) > 0, which is impossible, and a = —\TiB w \(r] + ()), which is incompatible 
with inequality (B.l). 
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Appendix C 


C.l Program Listings 

The Matlab routines related to the GUI comprise thousands of lines of code. For this reason, 
only essential routines that obtain thermodynamic properties, evaluate equilibrium points and 
find valve positions for a steady operating target are listed. 

C.1.1 Mixer Equilibrium Point 


mxr_eq.m 

“/Graphical illustration of the mixer operating point. See report. 

“/ Hanz Richter, PhD 

“/.NASA - Stennis Space Center, 2002 

“/Set pressure range for plotting: 

‘/.Make sure: 

“/.P>Ps 

“/.P1>P 

*/.Pg>P 

P= [5800: 100: 8000] ; 

for i=l :length(P) , 
if Pg>=2*P(i) 
f=2.423e-2*Pg; 
else 

f =2 . 857e-2*sqrt (Pg“2-P (i) “2) ; 
end; 

“/Mixer density from steady mass balance: 

mdot_l=l . 76e-2*Cvl*sqrt ( (Pl-P(i) ) *rho_l) ; 
mdot_g=f *Cvg*sqrt (Tg+460)*rho_g/Pg; 

rho_cv(i)=((mdot_l+mdot_g)/ (1 . 76e-2*Cve) ) “2/(P(i)-Ps) ; 

“/.Enthalpy of above locus : 

"/.Find pressure column: 
pcol=2*(l+(P(i)-2000)/100) ; 

"/.Find density range: 
minrho=hydro(61 ,pcol-l) ; 
maxrho=hydro ( 1 ,pcol-l) ; 

“/.Interpolate enthalpy 

if (rho_cv(i)>=minrho) & (rho_cv(i)<=maxrho) , 

h_th(i)=interpl (hydro( : ,pcol-l) ,hydro( : ,pcol) ,rho_cv(i) ) ; 

else 

h_th(i)=NaN ; 
end; 

“/.Finally, steady-state enthalpy: 
h_ss(i)=(mdot_l*hl+mdot_g*hg)/ (mdot_l+mdot_g) ; 

“/.Create surface by extrusion of the above: 

Z(: ,i)=h_ss(i)*ones(length(P) ,1) ; 

end; 


C.l. 2 Valve Positions for Prescribed Operating Condition 

“Zmxr_cv.m 
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‘/Calculates valve positions to achieve given exit flow and temperature, with 
‘/given mixer operating pressure. 

%****hydro_h and hydro_rho must be globalized and loaded in the workspace!!!***** 

°/,Hanz Richter, PhD 

'/.NASA - Stennis Space Center, 2002 

function [Cvl , Cvg , Cve] =mxr_cv (T , we , P) 

global hydro_h 
global hydro_rho 
global Ps 
global V 
global PI 
global Pg 
global hi 
global hg 
global Tg 
global T1 
global Cl 
global C3 
global rho_g 
global rho_l 


/ T: Desired temperature at exit valve outlet, deg. Fahrenheit 
‘/,we: Desired output mas sf low, lbm/sec. 

‘/,P: Desired mixer operating pressure (must ensure positive flow based on input and output pressures) 


'/, Check pressures for positive flow: 

if (P<=Ps) | (P>Pg) I (P>P1), 

error ( ’Reverse flow. Change pressures’); 
end; 

‘/.Find enthalpy from Ps and T by 2D interpolation on hydro_h (see hydro_desc.txt) 
if (Ps<=2000) | (Ps>=13500) | (T<=-400) | (T>=200) , 
error ( ’ Of f -range ’ ) ; 
end; 

Pvect= [2000 : 100 : 13500] ; 

Tvect= [-400 : 10 : 200] ; 

h=interp2(Tvect ,Pvect ,hydro_h’ ,T,Ps) 

‘/.Enthalpy accross exit valve is constant. Use it to find mixer density ( =at exit valve inlet) 

‘/.Generate interpolated enthalpy and density columns for given mixer pressure. 

‘/.Find bracketing pressures: 
plow_index=max(f ind(Pvect<=P) ) ; 
plow=2000+ (plow_index- 1 ) * 100 ; 
phigh=plow+100 ; 

hsegment=hydro_h( : , [plow_index,plow_index+l] ) ; 
rhosegment=hydro_rho( : , [plow_index,plow_index+l] ) ; 

‘/.Interpolated columns : 

hdata=interp2(Tvect , [plow.phigh] ,hsegment’ ,Tvect,P) ; 
rhodata=interp2(Tvect , [plow,phigh] ,rhosegment ’ ,Tvect,P) ; 

‘/.Find the density for h: 
if (h>max(hdata) ) I (h<min(hdata) ) , 
error ( ’ Of f -range ’ ) ; 
end; 

rho=interpl (hdata,rhodata,h) 

‘/.Find equilibrium flows (see report) 
wl=we*(h-hg)/ (hl-hg) ; 
wg=we*(hl-h)/ (hl-hg) ; 

‘/.Find exit Cv: 

Cve=we/ (C3*sqrt ( (P-Ps) *rho) ) ; 

‘/.Find gas Cv: 
if Pg>=2*P 

f=2.423e-2*Pg; 

else 

f =2 . 857e-2*sqrt (Pg“2-P“2) ; 
end; 

Cvg=wg/ (f *sqrt (Tg+460) *rho_g/Pg) ; 

‘/.Find liquid Cv: 

Cvl=wl/ (Cl*sqrt ( (Pl-P) *rho_l) ) ; 
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C.1.3 Thermodynamic Properties 

thermo_props . m 

°/,Hanz Richter, NASA 2002 


‘/.This function returns the temperature and pressure of nitrogen or hydrogen given density in lbm/ft“3 and 
'/, internal energy in Btu/lbm. 

'/, Final version (used to be called gettemp, getpress, and getpropsl,2,3, and 4. 

‘/.July 18, 2002 

function vect=thermo_props (rho , u , fluid) 


switch lower (fluid) 

case { ’parahydrogen’ , ’l’}, 

global hydro_u_data 
global hydro_p_data 
global hydro_t_data 
global hydro_master_range 

u_data = hydro_u_data; 

p_data = hydro_p_data; 

t_data = hydro_t_data; 

master_range = hydro_master_range ; 

rho_lower = .535; 

rho_step = 0.025; 

rho_upper = 6.31; 

u_lower = -88.8637; 

u.upper = 1559.8; 

nrows = 116; 

dens_interp_treshold=0 . 002 ; 
energy_interp_treshold=0 . 1 ; 

dens_perturb_step=0.02; '/these two control the error in original density 
dens_tol=0.25; 

case {’nitrogen’, ’2’}, 

global nitro_u_data 
global nitro_p_data 
global nitro_t_data 
global nitro_master_range 

u_data = nitro_u_data; 

p_data = nitro_p_data; 

t_data = nitro_t_data; 

master_range = nitro_master_range ; 

rho_lower = .05; 

rho_step = 0.5; 

rho_upper = 54.55; 

u.upper = 224 . 1 ; 

u_lower = -59.4982; 

nrows = 250; 

dens_interp_treshold=0 . 01 ; 
energy_interp_treshold=0 . 1 ; 

dens_perturb_step=0.02; '/these two control the error in original density 
dens_tol=0.25; 

otherwise , 

error (’ Invalid fluid’); 

end; 

if (rho>rho_upper) I (rho<rho_lower) , 

"/of f -range 

error ( ’Density off-range’) 
end; 

if (u>u_upper) I (u<u_lower) , 

'/of f -range 

error ( ’Energy off-range’) 
end; 


last_resort=0 ; 
interpolate_dens=l ; 
found=0; 
dens_error=0; 

while "found & dens_error<dens_tol, 

colnumber=floor((rho-rho_lower)/rho_step)+l; '/find left bracket column 

‘/.First priority: check if the density is "close" to one of the nominal bracketing densities 
if abs (rho- (rho_lower+ (colnumber-1) *rho_step) ) <=dens_interp_treshold , 
interpolate_dens=0 ; 
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elseif abs(rho-(rho_lower+(colnumber)*rho_step))<=dens_interp_treshold, ‘/check the next column as well 
interpolate_dens=0 ; 

colnumber=colnumber+l ; ‘/change the column 

end; 

if ~interpolate_dens , 

u_data_col=u_data( : .colnumber) ; 
p_data_col=p_data( : .colnumber) ; 
t_data_col=t_data( : .colnumber) ; 

‘/check if the density column contains a "close" energy match 
search_scope=master_range (colnumber) ; 

rownumber=nrows+max (f ind(abs (u_data_col (nrows-search_scope+l : nrows) -u) <=energy_interp_treshold) ) -search_scope ; 
if isempty(rownumber) I rownumber==nrows , ‘/.Now try for bracketing energy in the col: 

rownumber=nrows+max (f ind (u_data_col (nrows-search_scope+l : nrows) <=u) ) -search_scope ; 
if isempty(rownumber) I rownumber==nrows , 

interpolate_dens=l ; ’/See if interpolation helps it. 

else 

p=interpl(u_data_col(rownumber:rownumber+l) ,p_data_col(rownumber :rownumber+l) ,u) ; 
t=interpl(u_data_col(rownumber:rownumber+l) ,t_data_col (rownumber :rownumber+l) ,u) ; 
f ound=l ; 
end; 
else 

p=p_data_col(rownumber) ; 
t=t_data_col(rownumber) ; 
f ound=l ; 
end; 
end; 

‘/.Second priority: check if the interpolated density columns contain a "close" energy match 
if interpolate_dens, ’/.then colnumber should still be the left bracket 

colnumber=floor((rho-rho_lower)/rho_step)+l ; ‘/.find left bracket column 
f raction= (rho- (rho_lower+( colnumber- 1) *rho_step) ) /rho_step ; 

u_data_col=u_data( : , colnumber) +fract ion* (u_data( : , colnumber+l)-u_data( : .colnumber)) ; 
p_data_col=p_data( : , colnumber) +fract ion* (p_data( : , colnumber+l)-p_data( : .colnumber)) ; 
t_data_col=t_data( : , colnumber) +fract ion* (t_data( : ,colnumber+l)-t_data( : .colnumber)) ; 

‘/don’t use the zeros in the interpolation: 

search_scope=min(master_range (colnumber) ,master_range (colnumber+1) ) ; 

'/.Try locate the energy in u_data_col within search_scope 
'/.Remember data is stored at the bottom, and in ascending order 

'/ 0 _master_range tells how many rows are useful from the bottom up. 

rownumber=nrows+max (f ind (abs (u_data_col (nrows-search_scope+l : nrows) -u) <=energy_interp_treshold) ) -search_scope ; 
if isempty (rownumber) I rownumber==nrows , '/.Now try for bracketing energy in the interpolated col: 
rownumber=nrows+max (find (u_data_col (nrows-search_scope+l : nrows) <=u) ) -search_scope ; 
if isempty (rownumber) I rownumber==nrows , '/.In this case, go to last-resort procedure 

last_resort=l ; 
else 

p=interpl(u_data_col (rownumber :rownumber+l) , p_data_col (rownumber :rownumber+l) ,u) ; 
t=interpl(u_data_col (rownumber :rownumber+l) , t_data_col (rownumber :rownumber+l) ,u) ; 
f ound=l ; 
end; 
else 

p=p_data_col (rownumber) ; 
t=t_data_col (rownumber) ; 
f ound=l ; 
end; 
end; 

if last_resort, ‘/.As a last resort, change the density and start over 
rho=rho+dens_perturb_step ; 
dens_error=dens_perturb_step+dens_error ; 
end; 
end; 

if found, 

vect= [t p dens_error] ; 
else 
rho 
u 

dens_error 

error(’Pair really off-range’); 
end; 


C.2 Feedback Linearization Controller 

function w=f eedback_lin_control (x) 

*/,x contains [eve y,yd,ydd,P,T,pars,en] 

‘/,y and yd are arranged as [rho, u, exit_massf low] “T 
*/,P is pressure measurement 
‘/,T is temperature measurement 

'/.pairs is a vector containing [pl,hl,rhol,pg,hg,rhog,ps,Tg] 

‘/,en is an enable flag to avoid doing the calculations when running open-loop or with other 
'/.controllers, en must match this controller’s value for the calculations to be performed, 
y** ************************** ****** 

‘/.Feedback linearization controller flag = 2 
%* ************************** ******* 
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’/Therefore x=[ eve rho,u,we,rho_d,u_d,we_d,rho_d_d,u_d_d,we_d_d,P, T, pi, hi, rhol, 

*/. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

°/. pg i hg , rhog , ps , Tg] . 

’/. 16 17 18 19 20 

’/.The virtual control w is arranged as 

'/,w= [Cvl ; Cvg ; v] , where v is integrated to give eve. 


global V Clp C2p C3p C3pp C4p; 
if x(21) == 2 


cve=x(l) ; 

'/.Calculate E 

auxl=Clp*sign(x(13)-x(ll) )*sqrt (abs(x(13)-x(ll) )*x(15)) ; ’/.this is fl 
aux2=sign(x(ll)-x(19) )*sqrt (abs(x(ll)-x(19) )*x(2)) ; 
aux3=x ( 2 ) *ge tpart i alrho (x(ll),x(12))+x(ll)-x(19); 

if x(16)<2*x(ll) , 

aux4=C4p*sign(x(16)-x(ll) )*sqrt (abs(x(16) ~2-x(ll) “2) ) ; 
else 

aux4=C2p ; 
end; 


E=[auxl aux4 0;auxl*(x(14)-x(3))/x(2) aux4*(x(17)-x(3))/x(2) 0; -aux3*V*cve*C3p*auxl/ (2*aux2) -aux3*V*cve*aux4*C3p/(2*aux2) -V*C3p*aux2] 
'/.Calculate D 

D= [C3p*cve*aux2;C3pp*cve*aux2*x(ll)/x(2) “2; -V*cve“2*C3p“2*aux3/2] ; 

’/.Control gains 
gam=diag( [10 10 5]); 

’/.Calculate the virtual control input 
w=inv(E)*(x(8: 10)-gam*(x(2:4)-x(5:7))-D) ; 
else 

w= [0 0 0]; 

end 


C.3 Sliding Mode Controller 

‘/.smctrl .m 

’/.This function calculates the sliding mode controller 
’/.Hanz Richter, NASA 2002 


function sliding_control = smctrl (x) 


global V c Ps rhol rho2 PI P2 R hi h2 Cv k taul tau2 tau3 betal beta2 beta3 


f l=x(l) ; 
f2=x(2) ; 
fe=x(3) ; 
rho=x(4) ; 
P=x(5) ; 
Cvl=x(6) ; 
Cv2=x(7) ; 
Cve=x(8) ; 
rho_d=x(9) ; 
rho_dd=x(10) ; 
rho_ddd=x(ll) 
P_d=x(12) ; 
P_dd=x(13) ; 
P_ddd=x(14) ; 
w_d=x(15) ; 
w_dd=x(16) ; 


’/.rho desired 
’/.rho desired derivative 
'/.rho desired second derivative 
°/,P desired 

°/,P desired derivative 

°/,P desired second derivative 

'/.flow desired 

'/.flow desired derivative 


krho=100; ’/.sliding manifold coefficients 

kP=100; 

kw=100; 

etal=120; ’/.sliding gains 

eta2=2400e3; 
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eta3=120; 


phil=10; 

phi2=le4; 

phi3=10; 


f=[fl f 2 -fe] ’ ; 

Cvect= [Cvl ; Cv2 ; Cve] ; 

Df lDp=-0 . 5*c*rhol/sqrt (abs (Pl-P) *rhol) ; 

Df 2Dp=-0 . 5*c*rho2/sqrt (abs (P2-P) *rho2) ; 

Df eDp=0 . 5*c*sqrt (rho) / sqrt (abs (P-Ps) ) ; 

Df Drho= [0 ; 0 ; -0 . 5*c*sqrt (abs (P-Ps) /rho) *sign(P-Ps) ] ; 

Df Dp= [Df lDp ; Df 2Dp ; -Df eDp] ; 

Ap=diag( [R*hl/Cv R*h2/Cv k*P/rho]); 

Ac=diag( [-1/taul -l/tau2 -l/tau3]); 

Bc=diag( [betal beta2 beta3] ) ; 

'/.calculation of sliding functions 
sl=rho_dd-f ’ *Cvect/V+krho* (rho_d-rho) ; 
s2=P_dd-f ’ *Ap*Cvect/V+kP* (P_d-P) ; 
s3=kw*(w_d-Cve*f e) ; 

'/.Flow Derivatives 

Df Dt= ( (Cvect ’ *Ap*f ) *Df Dp+ (Cvect ’ *f ) *Df Drho) /V ; 

‘/.Matrix derivative 
DApDt=zeros(3,3) ; 

DApDt (3,3) =k*Cvect ’ * (Ap-P*eye (3) / rho) *f /rho/V ; 

'/.Sliding Controller Calculations 

Gammal=rho_ddd+krho*rho_dd+etal*sat (s 1/phi 1) -Cvect ’* (DfDt+ (Ac+krho*eye (3,3) )*f )/V ; 
Gamma2=P_ddd+kP*P_dd+eta2*sat (s2/phi2) -Cvect ’ * (Ap*Df Dt+ (Ap* (Ac+kP*eye (3,3)) +DApDt ) *f ) /V ; 
Gamma3=kw* (w_dd+Cve*f e/tau3+Cve*Df Dt (3) ) +eta3*sat (s3/phi3) ; 

G=zeros(3,3) ; 

G (3 , 3) =-kw*bet a3 ; 

syst_matrix= [f ’ *Bc/V ; f ’ *Ap*Bc/V ; f ’ *G] ; 
epsilons=inv (syst_matrix) * [Gammal ; Gamma2 ; Gamma3] ; 
sliding_control = [epsilons ; si ; s2 ; s3] ; 
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